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Tornado  -  Borne  Missile  Speeds 


Emil  Simiu  and  Martin  Cordes 

•■      At  the  request  of  the  United  States  Nuclear  Regulatory  Connnission  (NRC)  the  National 
Bureau  of  Standards  (NBS)  has  carried  out  an  independent  investigation  into  the  question 
of  tornado-borne  missile  speeds,  with  a  view  to  assisting  NRC  in  identifying  pertinent 
areas  of  uncertainty  and  in  estimating  credible  tornado-borne  missile  speeds  -  within  the 
limitations  inherent  in  the  present  state  of  the  art.     The  investigation  consists  of  two 
parts:     (1)  a  study,  covered  in  this  report,  in  which  a  rational  model  for  the  missile 
motion  is  proposed,  and  numerical  experiments  are  carried  out  corresponding  to  various 
assumptions  on  the  initial  conditions  of  the  missile  motion,  the  structure  of  the  tornado 
flow,  and  the  aerodynamic  properties  of  the  missile;   (2)  a  theoretical  and  experimental 
study  of  tornado-borne  missile  aerodynamics,  conducted  by  Colorado  State  University  (CSU) 
under  contract  with  NBS,  to  be  covered  in  a  separate  report  by  CSU.     In  the  present 
report,  the  factors  affecting  missile  motion,  and  their  influence  upon  such  motion,  are 
examined.     Information  is  provided  on  a  computer  program  developed  for  calculating  missile 
speeds.     Maximum  speeds  for  a  number  of  specified  potential  tornado-borne  missiles  are 
presented,  corresponding  to  a  set  of  assumptions  believed  by  the  writers  to  be  reasonable 
for  design  purposes.     It  is  pointed  out  that  higher  speeds  are  conceivable  if  it  is 
assumed  that  certain  circumstances,  examined  in  the  body  of  the  report,  will  obtain.  It 
is  the  judgment  of  the  writers  that  the  probabilities  of  occurrence  of  such  higher  speeds 
for  any  given  tornado  strike  are  low.     More  than  qualitative  estimates  of  such  probabilities 
are,  however,  beyond  the  scope  of  this  investigation. 
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1.  INTRODUCTION 


To  ensure  the  safety  of  nuclear  power  plants  in  the  event  of  a  tornado  strike  it  is 
required  that,  in  addition  to  the  direct  action  of  the  wind  and  of  the  moving   ambient  pressure 
field,  the  designer  consider  the  impart  of  tornado-borne  missiles,  i.e.,  of  objects 
moving  under  the  action  of  aerodynamic  forces  induced  by  the  tornado  wind.     It  is,  there- 
fore, necessary  that  estimates  be  made  of  the  speeds  attained  by  potential  missiles  under 
tornado  wind  conditions  specified  for  the  design  of  nuclear  power  plants. 

Several  such  estimates  have  been  reported  so  far   [1,  2,  3,  4,  5,  6,  7].     In  certain 
instances  the  differences  between  various  estimates  are  very  large;   for  example,  the 
predicted  speeds  of  a  utility  pole  [13-in  (0.33  m)  diameter  and  a  35-ft  (10.7  m)  length] 
predicted  by  Refs.  5,  6,  1,  2,   4  and  3  are  16.5  m/s,  30.5  m/s,  41.2  m/s,  42.7  m/s,  52.2 
m/s  and  54.9  m/s,  respectively.     For  any  given  missile,  the  kinetic  energy  associated 
with  translational  motion  is  proportional  to  the  square  of  the  speed;  in  the  case  of  the 
utility  pole,  therefore,  the  ratio  of  the  largest  to  the  smallest  of  the  kinetic  energy 
estimates  of  Ref.  1  through  Ref.   6  is  greater  than  10.     For  other  missiles,  the  ratio 
varies  from  almost  2  in  the  case  of  a  1-in  (2.54  cm)  diameter  steel  rod  to  over  5  in  the 
case  of  a  4000  lb   (18,000N)  standard  automobile. 

In  most  cases,  such  large  discrepancies  are  the  consequence  of  differences  between 
basic  assumptions  used  in  the  various  estimation  procedures.     These  assumptions  include: 

-  the  initial  conditions  of  the  problem,  i.e.,  the  initial  position  of  the  object 
with  respect  to  the  ground  and  to  the  tornado  center,  and  the  initial  velocity 
of  the  object. 

the  detailed  features  of  the  wind  flow  field. 

-  the  aerodynamic  characteristics  of  the  object  (which,  in  most  cases  of  practical 
interest,  is  a  bluff  body) . 

Differences  between  the  various  sets  of  basic  assumptions  used  in  estimating  tornado- 
borne  missile  velocities  may  be  ascribed,  in  part,  to  the  probabilistic  nature  of  the 
problem.     Indeed,  for  any  given  tornado  wind  field,   the  initial  conditions  constitute  a 
set  or  random  variables  with  a  very  large  number  of  possible  values,  the  choice  of  which 
is  not  unique.    More  importantly,  however,  such  differences  are  a  consequence  of  serious 
uncertainties  regarding  both  the  structure  of  the  tornado  flow  and  the  aerodynamic  behavior 
of  the  potential  missile. 
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At  the  request  of  the  United  States  Nuclear  Regulatory  Commission  (NRC) ,  the  National 
Bureau  of  Standards  has  carried  out  an  independent  Investigation  into  the  question  of 
tornado-borne  missile  speeds,  with  a  view  to  assisting  NRC  in  identifying  pertinent  areas 
of  uncertainty  and  in  estimating  credible  tornado-borne  missile  speeds — within  the  limitations 
inherent  in  the  current  state  of  the  art. 

The  objectives  of  the  investigation  were  defined  as  follows: 

1.  To  select  a  rational  model  for  the  tornado-borne  missile  motion. 

2.  To  develop  a  computer  program  based  on  this  model,  capable  of  computing  missile 
trajectories  and  velocities  for  any  specified  initial  conditions,  tornado  wind 
speed  model,  and  assumptions  regarding  the  drag  force  acting  on  the  missile. 

3.  To  calculate,  for  a  specified  set  of  initial  conditions  and  for  a  specified 
wind  speed  model,  the  trajectories  and  velocities  corresponding  to  a  number  of 
specified  potential  missiles. 

4.  To  determine,  in  a  number  of  representative  cases,  the  sensitivity  of  the 
calculated  results  to  changes  in  the  assumed  initial  conditions  or  in  the 
assumed  tornado  wind  speed  model. 

5.  To  obtain  and  interpret  information,  based  on  wind  tunnel  tests,  regarding  the 
aerodynamic  behavior  of  various  potential  missiles,  i.e.,  drag  coefficients  for 
various  missile  motions,  including  motion  in  the  tumbling  mode. 

6.  From  the  results  of  items  3,  4  and  5,  suggest  credible  speeds  of  selected 
potential  tornado-borne  missiles,  compatible  with  the  current  state  of  the  art 
in  nuclear  power  plant  design. 


2.     MODEL  FOR  THE  TORNADO-BORNE  MISSILE  MOTION 


The  motion  of  an  object  may  be  described  in  general  by  solving  a  system  consisting 
of  three  equations  of  balance  of  momenta  and  three  equations  of  balance  of  moments  of 
momenta.     In  the  case  of  a  bluff  body,  one  major  difficulty  in  writing  these  six  equations 
is  that  the  aerodynamic  forcing  functions  are  not  known. 

It  is  possible  to  measure  in  the  wind  tunnel  aerodynamic  forces  and  moments  acting 
on  a  bluff  body  under  static  conditions  for  a  sufficient  number  of  positions  of  the  body 
with  respect  to  the  mean  direction  of  the  flow.    On  the  basis  of  such  measurements,  the 
dependence  of  the  forces  and  moments  on  position,  and  corresponding  aerodynamic  derivatives, 
can  be  obtained.     Aerodynamic  forces  and  moments  can  then  be  calculated  following  the 
well-known  pattern  used  in  airfoil  theory;  for  example,  if  an  airfoil  has  a  time-dependent 
vertical  motion  h(t)  in  a  uniform  flow  with  velocity  U,  and  if  the  angle  of  attack  is 
a  =  const,  the  lift  coefficient  is 

r    -  '^^L     ,    _^  1  dh, 

^L  -  d^  U    dF)  (1) 

This  procedure  for  calculating  aerodynamic  forces  and  moments  is  valid  if  the  quasi- 
steady  assumption  [Ref.  8,  p.  192]   is  acceptable  and  if  the  body  concerned  behaves 


aerodynamically  like  an  airfoil  -  i.e.,  if  the  body  is  streamlined  and  if  no  flow  separation 
occurs.     However,  in  the  case  of  unconstrained  bluff  bodies  moving  in  a  wind  flow  the 
validity  of  such  a  procedure  remains  to  be  demonstrated. 

In  the  absence  of  a  satisfactory  model  for  the  aerodynamic  description  of  the 
missile  as  a  rigid   (six -degrees-of -freedom)  body,  it  is  customary  to  resort  to  the  alternative 
of  describing  the  missile  as  a  material  point  acted  upon  by  a  drag  force 


D  =  1/2  p  C^A  r        -         I    (v^  -  ■  • 


where  p  =  air  density,  v    =  wind  velocity,  v    =  missile  velocity,  A  is  a  suitably  chosen 

w  .  M 

area  and  C^^  is  the  corresponding  drag  coefficient. 

This  model  is  reasonable  if,  during  its  motion,  the  missile  either  (a)  maintains  a 
constant  or  almost  constant  attitude  with  respect  to  the  relative  velocity  vector  v 

w 

v„»  or  (b)  has  a  tumbling  motion  such  that,  with  no  significant  errors,  some  mean  value 
M   

of  the  quantity  C^^A  can  be  used  in  the  expression  for  the  drag  D.     The  assumption  of  a 
constant  body  attitude  with  respect  to  the  flow  would  be  credible  if  the  aerodynamic 
force  were  applied  at  all  times  exactly  at  the  center  of  mass  of  the  body — which  is 
highly  unlikely  in  the  case  of  a  bluff  body  in  a  tornado  flow — ,  or  if  the  body  rotation 
induced  by  a  non-zero  aerodynamic  moment  with  respect  to  the  center  of  mass  were  inhibited 
by  aerodynamic  damping  forces  intrinsic  in  the  body-fluid  system.     The  question  thus 
arises  as  to  whether  such  stabilizing  forces  may  be  expected  to  be  present. 

It  is  of  interest  at  this  point  to  mention  certain  experimental  results — obtained  in 
studies  of  bridge  deck  aerodynamic  stability — which  provide  useful  insights  into  the 
question  at  hand.     Consider  a  body  restrained  by  four  springs  of  equal  stiffness,  immersed 
in  a  horizontal  flow  (Fig.  1) ,  and  subjected  to  an  impulse  which  produces  angular  oscillations 
6(t)  about  the  position  of  equilibrium  [9].     In  the  case  of  an  airfoil  with  a  sufficiently 
small  angle  of  attack  so  that  flow  separation  does  not  occur,  the  aerodynamic  damping, 
which  is  proportional  to  the  quantity  denoted  by  H*  in  Ref .  9,  is  positive.     This  implies 
that  the  flow  will  contribute  ,  along  with  the  viscous  damping  inherent  in  the  springs, 
to  the  damping  out  of  the  oscillations.     On  the  other  hand,  for  bluff  bodies,  at  high 
velocities  of  the  flow  and  for  vanishing  values  of  the  spring  stiffness,  the  aerodynamic 
damping  is  negative  in  the  large  majority  of  the  cases  tested  [9] 

These  results  suggest  that,  in  general,  no  stabilizing  effect  by  the  flow  can  be 
expected  to  inhibit  the  tumbling  of  bluff  bodies.     The  assumption  that  potential  tornado- 
borne  missiles  will  tumble  during  their  motion  appears  therefore  to  be  reasonable.  It 
will  be  this  assumption  that  will  fee  used  in  this  work. 
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Assuming  that  Eq.  2  is  valid  and  that  the  average  lift  force  vanishes  under  tumbling 
conditions,  the  motion  of  the  missile  viewed  as  a  one-degree-of-f reedom  system  is  governed 
by  the  relation: 


'^^M           ,  ,„      C  A     I  -  -     I     -       -           -  (3) 

-r—    =  -  1/2  p  _D_       V  -  V        (v^  -  V  )  -  gk 

dt                                   '     M  w  '       M  w 

m 

where  g  =  acceleration  of  gravity,  k  =  unit  vector  along  the  vertical  axis  and  m  =  mass  of 
missile. 

It  follows  from  Eq.  3  that  for  a  given  flow  field  and  initial  conditions,  the  motion 

depends  only  upon  the  value  of  the  parameter  C  A/m.     Throughout  this  work,  the  area.  A,  will 
2 

be  expressed  in  m    and  the  mass,  m,  in  kg.     To  transform  the  parameter  C  A/w,  where  A  is 
2 

expressed  in  ft    and  the  weight  in  lb,  into  the  parameter  C  A/ m,  where  A  is  expressed  in 
2 

m    and  the  mass,  m,  in  kg,  the  following  relations  are  used: 

-2-    =  1  If-  =  .205  f-  ^^^^ 

w  lb  m  kg 


C^A  2  C^A  ^2 

^    =l2_  ->  =  4.902  1^-  (4b) 

m  kg  w  lb 

3.     COMPUTER  PROGRAM  FOR  CALCULATING 
TORNADO-BORNE  MISSILE  TRAJECTORIES  AND 
VELOCITIES 

To  calculate  and  plot  trajectories  and  velocities  of  tornado-borne  missiles,  a 

computer  program  was  developed  in  which  the  assumed  models  for  the  tornado  wind  field  and 

the  drag  coefficients  are  specified  by  specialized  subroutines  (details  of  such  models 

are  given  in  subsequent  sections) . 

In  Eq.  3  the  components  of  the  missile  velocity  v,,  and  of  the  wind  velocity  v  must 

_M  w 

be  referred  to  an  absolute  frame.     The  wind  velocity  v^  is  usually  specified  as  a  sum  of 

two  parts.     The  first  part  represents  the  wind  velocity  of  a  stationary  tornado  vortex 

and  is  referred  to  a  cylindrical  system  of  coordinates.     The  second  part  represents  the 

translation  velocity  of  the  tornado  vortex — or,  equivalently ,  of  the  cylindrical  system 

of  coordinates — with  respect  to  an  absolute  frame  of  reference.     The  transformations 

required  to  represent  v    in  an  absolute  frame  are  derived  in  Appendix  A  and  are  incorporated 

w 

in  the  computer  program.     Documentation  on  the  computer  program  and  a  sample  input  and 
output  are  given  in  Appendix  C.     The  program,  which  is  written  in  ANSI  Fortran  language, 
may  be  obtained  on  tape  from  the  National  Technical  Information  Service,  Springfield, 
Virginia  22151. 

For  the  particular  case  of  a  parallel  flow,  an  analytical  solution  to  the  problem  of 
the  missile  trajectory  can  be  easily  obtained.     This  solution,  which  can  be  found  in 
Appendix  B,  has  been  used  to  test  the  program,  in  which  the  subroutine  describing  the 
wind  field  was  modified  to  represent  a  parallel  flow. 


4.     NUMERICAL  COMPUTATIONS 


It  was  previously  noted  that,  for  a  given  flow  field  and  for  given  initial  conditions, 
the  missile  motion  depends  only  upon  the  value  of  the  parameter  C^A/m.     In  this  section, 
numerical  results  will  be  presented  which  show  the  effect  of  this  parameter  on  the  maximum 
horizontal  missile  speed.    The  calculations  are  based  on  the  following  assumptions: 

1.  The  parameter  C^A/m  is  constant  during  the  missile  flight. 

2.  The  tornado  wind  field  may  be  represented  by  a  vortex  translating  with  a 

2 

uniform  velocity  v^  along  an  axis  denoted  by  Ox  (Fig.  2).     Let  v„,  v„ ,  v    and  v       =  (v„ 
2  R      9      z  rot  R 

+  V    )        denote  the  radial,  tangential,  vertical  and  horizontal  velocity  in  the  vortex 

6 

flow,  respectively.     It  is  assumed  that 

^rot  =|-v  0<r<R  (5) 

R        m  —      —  ra 

m 

R 

rot  =  —    V  R<r<<»  (6) 

r        m  m  —  — 

where  v    is  the  maximum  horizontal  velocity  in  the  vortex  flow,  r  is  the  radius  (distance 
m 

from  the  vortex  center)  and  R    is  the  value  of  r  at  which  this  velocity  is  attained. 

m 

Eqs.  5  and  6  are  similar  to  descriptions  of  the  flow  proposed  in  Ref .  10  (p.  135). 
Furthermore,  it  is  assumed  that        =  1/2  v^  and  v^  =  0.67  v^ ,  i.e.,  (see  Fig.  2). 

1 

V  =  —  V 

R  rot  (7) 

2 

^     ''e       /"s-'^rot  (8) 

V  =    V  (9) 

z  3/^ 

The  type  of  model  just  described  is  referred  to  as  the  Rankine  vortex  and  appears  to 

be  a  reasonable  representation  of  tornado  flows.     Estimates  based  on  field  observations 

suggest  that  it  is  reasonable  to  assume — as  is  done  in  this  model — that  R  is  independent 
of  height  [Ref.  10,  p.  131]. 
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The  following  values  for  the  tornado  wind  field  parameters  were  used  in  the  calculations: 


Table  1  -  Values  of  v„,  v    and  R    Used  in  Numerical  Calculations 
T      m  m 


v_ 

V 

V  = 

R 

i 

torn 

m  1 

TH 

Tornado 

m/s* 

mph 

m/s* 

mph 

m/  s* 

mph 

m* 

ft 

Type 

1 

31 

70 

130 

290 

161 

360 

46 

150 

2 

27 

60 

107 

240 

134 

300 

46 

150 

3 

22 

50 

85 

190 

107 

240 

46 

150 

4 

31 

70 

146 

325 

177 

395 

46 

150 

5 

27 

60 

120 

288 

147 

348 

46 

150 

6 

22 

50 

95 

213 

117 

263 

46 

150 

*  Approximately 

The  values  given  in  Table  1  for  tornado  types  1,  2,  3  are  suggested  in  Ref.  11  as  providing 

an  acceptably  low  level  of  failure  if  used  in  the  design  of  nuclear  power  plants.  The 

values  for  tornado  type  4,  5  and  6  were  included  for  the  purpose  of  studying  the  effect 

upon  missile  velocity  of  relatively  small  increments  in  the  value  of  v^. 

3.     The  assumed  initial  conditions  are:     x(0)  =  R  ,  y(0)  =  0,  z(0)  ,=  40m,  v„  (0)  =  0, 

m  Mx 

v.,  (0)  =  0,  v.,  (0)  =  0  at  time  t  =  0,  where  x,  y,  z  are  the  coordinates  of  the  missile 
My  Mz 

(i.e.,  of  its  center  of  mass)  and  v.,  ,  v^,  ,  v,,    are  the  missile  velocity  components  alone 

Mx      My      Mz  J        r  b 

the  X,  y,  z  axes  (Fig.  2),    Also,  at  time  t  =  0  the  center  of  the  tornado  vortex  coincides 

with  the  origin  0  of  the  coordinate  axes.     The  effect  of  assuming  initial  conditions  different 

from  those  indicated  is  examined  in  the  next  section  of  this  report. 

The  dependence  upon  the  parameter  C^A/m  of  the  maximum  horizontal  missile  speed  calculated 

in  accordance  with  assumptions  1  through  3  is  represented  in  Fig.  3  for  tornado  types  1 

through  6  as  defined  in  Table  1.     In  Fig.  3  v          =  v^  +  v      (Table  1). 

torn        i  m 

5.     SENSITIVITY  OF  CALCULATED  RESULTS 

TO  CHANGES  IN  THE  ASSUMED  INITIAL 
CONDITIONS  OR  IN  THE  ASSUMED  TORNADO 
WIND  SPEED  MODEL 

5.1    Changes  in  the  Assumed  Initial  Position  of  the  Missile. 

For  flows  with  v    =  146  m/s,  R    =  46  m,  and  v„  =  31  m/s,  the  maximum  horizontal 
mm  1 

speeds  of  missiles  with  C^A/m  =  0.001  and  C^A/m  =  0.01  were  calculated  using  the  initial 
positions  shown  in  Table  2.     Except  for  the  initial  positions,  the  results  of  Table  2  are 
based  on  the  assumptions  described  in  the  preceding  section.     It  is  seen  from  Table  2  that 
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Table  2  -  Maximum  Horizontal  Missile  Speeds,  v^^^  ,  (m/s) 
Corresponding  to  Various  Initial  Positions 


Initial 

Position 

Cj^A/m 

(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

x(0) 

y(0) 

(J  .  UUl 

U .  U.L 

(meters) 

(a) 

o- 

46 

0 

10 

65 

(b) 

G- 

23 

0 

18 

93 

(c) 

69 

0 

9 

48 

(d) 

o- 

-46 

0 

18 

82 

(e) 

0 

46 

16 

68 

(f) 

0 

23 

20 

84 

(g) 

0 

-23 

35 

50 

(h) 

0 

-46 

54 

70 

1.  Arrows  in  column  (2)  represent  direction  of  translation  velocity 

2.  Assumed  elevation  of  missile  at  time  t  =  0:     z(0)  =  40m  in  all  cases. 


Table  3  -  Maximum  Horizontal  Missile  Speeds,  (m/s)  Corresponding  to 

Initial  Elevations  2(0)  =  10m. 
z(0)  =  20m  and  z(0)    =     40m  (C  A/m  =  0.001) 


Initial 

Position 

z  (0) 

(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

^(0) 

y(0) 

10m 

20m 

40m 

(meters) 

(a) 

o 

-> 

46 

0 

6 

7 

10 

(b) 

G 

-> 

23 

0 

2 

8 

18 

(c) 

-> 

69 

0 

9 

9 

9 

(d) 

-> 

-46 

0 

14 

16 

18 

(e) 

0 

-»■ 

0 

46 

3 

9 

16 

(f) 

-»■ 

0 

23 

12 

16 

20 

(g) 

Q 

-> 

0 

-23 

19 

29 

35 

(h) 

Q 

->■ 

0 

-46 

33 

42 

54 
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the  initial  position  used  in  the  calculations  of  the  preceding  section,   (position  (a)  of 
Table  2)  does  not  result  in  the  largest  possible  maximum  horizontal  missile  speeds. 
It  is  also  noted  that  the  initial  positions  to  which  there  correspond  the  largest  possible 
missile  speeds  depend  upon  the  value  of  C^A/m.     For  example,  for  C^A/m  =  0.01,  that  initial 
position  is  (b ) ;  for  C^k/m  =  0.001,  it  is  (h)   [see  Table  2]. 

As  indicated  in  Note  2  of  Table  2,  the  initial  elevation  assumed  in  the  calculations 
was  z(0)  =  40m.     If  the  weight  of  the  missile  is  smaller  than  the  upward  drag  induced 
by  the  vertical  wind  velocity  component,  then  the  calculated  missile  velocities  at  time  t 
are  independent  of  2(0) .     (This  is  the  case  because,  in  the  assumed  tornado  model,  the 
flow  field  is  invariant  with  z) .    However,  if  the  missile  weight  exceeds  the  upward  drag, 
i.e.,  if  the  missile  moves  downward, — as  in  the  case  of  column  (5)  of  Table  2 —  the 
interval  between  time  t  =  0  and  the  time  at  which  the  missile  hits  the  ground  decreases  as 
z(0)  decreases.    Therefore,  the  maximum  missile  speed  may  decrease  if  lower  values  of  the 
initial  elevation  z(0)  are  assumed.     Table  3  lists  speeds  calculated  using  the  assumptions 
z(0)  =  10m  and  z(0)  =  20m,  all  other  values  of  the  various  parameters  being  the  same  as 
for  column  (5)  of  Table  2.     For  comparison  the  calculated  speeds  in  column  (5)  of  Table  2 
were  also  included  in  Table  3. 

Missile  speeds  were  also  calculated  corresponding  to  various  initial  elevations 
0  <  z  (0)  <  10m.     It  is  of  interest  to  note  that  the  maximum  horizontal  speed  of  a  missile 
with  CjjA/m  =  0.001  starting  the  motion  from  position  (h)  is  relatively  high  even  for  low 
values  of  the  initial  elevation  z(0).     For  example,  for  z(0)  =  3m  and  z(0)  =  5m,  the  maximum 
missile  speeds  are  23  m/s  and  27  m/s,  respectively. 

Calculations  were  also  carried  out  for  the  horizontal  distances  traveled  by  the 
missiles.    For  example,  for  C^A/M  =  0.001  and  the  initial  position  (h) ,  the  horizontal 
distances  corresponding  to  z(0)  =  3m,  5m,  10m,  20m,  and  40m  are,  approximately,  20m,  30m, 
50m,  90m  and  160m,  respectively. 

5.2    Changes  in  the  Assumed  Initial  Velocity  of  the  Missile. 

The  results  given  in  the  preceding  sections  are  based  on  the  assumption  that  the 

initial  velocity  of  the  missiles  is  zero.     If  the  missile  is  injected  in  the  flow,  for 

example  by  an  explosion,  this  assumption  is  no  longer  valid.     However,  all  other  conditions 

being  equal,  a  non-zero  intial  velocity  does  not  necessarily  result  in  a  maximum  missile 

velocity  higher  than  that  corresponding  to  zero  initial  velocity.     This  is  illustrated 

in  Table  4,  in  which  type  1  tornado  (see  Table  1),  and  the   conditions  ^^^^^(0)  ^  0, 

Vw     (0)  =  0,  v.,     (0)  =  0  were  assumed. 
My  Mz 

Table  4  -  Maximum  Horizontal  Missile  Speeds,  v™^^,   (m/s)  Corresponding  to  Various  Initial  Velocities 


C^k/m  =  0.001 

C^A/m  =0.01 

(2) 

(3) 

(4) 

(5) 

(7) 

(8) 

(9) 

(10) 

x(0) 

y(0) 

(meters /second) 

(meters) 

0 

10 

20 

0 

10 

20 

(a) 

46 

0 

8 

9 

20 

62 

58 

53 

(g) 

0 

-23 

35 

45 

35 

63 

59 

59 
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5.3  Changes  in  the  Assumed  Tornado  Wind  Speed  Model. 


5.3.1  Tornado  Described  by  Eqs.  5  through  9  with  v    =  0. 


Table  5  lists  maximum  horizontal  missile  speeds  for  the  same  parameter  values  and 
initial  conditions  as  those  used  in  Table  2,  except  that  the  translation  velocity  of  the 
tornado  vortex  is  zero,  rather  than  31  m/s.     It  is  seen  by  comparing  Tables  2  and  5  that 
the  calculated  speeds  are,  on  the  average,  higher  in  the  case  v^  =  31  m/s,  although  for 
some  initial  positions  the  reverse  is  true. 


Table  5  -  Maximum  Horizontal  Missile  Speeds,  v^ 
Corresponding  to  Various  Initial  Positions  (v    =  0) 


Initial  Position 

Cj^A/m 

(1) 

(2) 

x(0) 

y(0) 

(2) 

(3) 

(meters) 

0.001 

0.01 

(a) 

o 

46 

0 

30 

65 

(b) 

23 

0 

7 

61 

(c) 

8- 

69 

0 

19 

60 

(d) 

-46 

0 

30 

65 

(e) 

6 

0 

46 

30 

65 

(f) 

0 

23 

8 

61 

(g) 

0 

-23 

8 

61 

(h) 

0 

-46 

30 

65 

5.3.2  Tornado  Vortex  Models  Based  on  Data  Obtained  During  the  Dallas  Tornado  of 
April  2.  1957  [13]. 

In  this  family  of  models,  the  velocities  in  the  tornado  vortex  depend  upon  the  quantities 

=  radius  of  maximum  tangential  velocity  at  elevation  z  =  0,  Vg^  =  maximum  tangential 

velocity  and  R*(z)  =  radius  beyond  which  the  radial  and  vertical  velocity  components 

vanish.     In  addition,  the  flow  depends  on  three  parameters  k^^,  k^,  k^,  as  will  be  shown 

subsequently.    The  radius  of  maximum  tangential  velocity  at  elevation  z  is  denoted  by 

R    (z)  and  is  assumed  to  be 
m 

R    (z)  =  R    +  k,  z  0  <  z  <  60m  (10) 

m  m       1  ^  ' 


R     (z)  =  R    +  60  k,  z  >  60m  (11) 

m  m  1  — 

It  is  seen  that  the  parameter  k^^  is  a  measure  of  the  extent  to  which  the  tornado 
vortex  deviates  from  a  cylindrical  shape.     If  k^^  =  0  the  radius  of  maximum  tangential 
velocity  is  invariant  with  height,  as  in  the  model  described  by  Eqs.  5  through  9.  An 
assumption  commonly  employed  in  tornado-borne  missile  speed  investigations  is  k^  =  0.45 
[see,  for  example,  Refs.  3  and  12],     It  is  further  assumed  [3,  12]: 
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R*  (z)  =  CR^  (z)  (12) 

where  C  Is  a  coefficient  which  depends  upon  Vg^.     In  this  work,  it  was  assumed  C  =  2.35 

for  V.    =  130  m/s,  C  =  2.10  for  v„    =  107  m/s  and  C  =  1.80  for  v„    =  85  m/s.     These  values 
om  t)m  bm 

are  similar  to  those  used  in  Ref.  3. 

Let  r  denote  the  distance  to  the  center  of  the  tornado  vortex.    The  components 
of  the  velocity  in  the  tornado  vortex  are  assumed  to  be 


0  <  z  <  60m 


V    -      k     R*  (z)  -  r  „^  ,  , 

■R  "        2    R*         -  R    (2)  r         r  <  R*  (2)  (13) 
m 


"  0  r  >  R*(z)  (14;, 


\  =  r-fir  ^m     •         ^  <  ^^^^ 


m 


—  m 


R 


V.  r  >  R    (z)  (16) 


r  9m 


m 


\  '  S  R*  [I]  :  R    (z)  ^     r  r<K*(z)  (17) 


«=  0  r  >  R*  (z)  (18) 


60  m  <  z  <  240m 


,     R*  (z)  -  r  240  -  z  ,  ^  /-.ox 

^  R*  (z)  -  R^  (z)  ^  '  ^  ^* 


=  0  r  >  R*  (z)  (20) 


v_  = 


Tl^  ^em  ^  (21) 

\  (z) 

^m  >  \  (22) 


-z  "  -  lfo>   (^3  R*  (z!  :  R  (z)  ^°      r>     ^  <  ^*  (23) 


=  0  r  >  R*  (z)  (24) 
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Table  6  -  Maximum  Horizontal  Missile  Speeds  v"^"^.   (m/s)  Corresponding  to  Various 
Parameter  Values 


Case 

^em 

^2 

"3 

Cj^A/m 

max 
\ 

(m/s) 

(m/s) 

m-^/kg 

(m/s) 

1 

85 

22 

.45 

1 

1 

0.025 

90* 

2 

85 

22 

.45 

1 

0 

0.025 

83* 

3 

85 

22 

0 

1 

0 

0.025 

77 

4 

85 

22 

0 

1 

1 

0.025 

79 

5 

130 

31 

.45 

1 

1 

0.025 

90* 

6 

130 

31 

.45 

1 

0 

0.025 

126* 

7 

130 

31 

0 

1 

0 

0.025 

80* 

8 

130 

31 

0 

1 

1 

0.025 

75* 

9 

130 

31 

.45 

1 

1 

0.01 

91* 

10 

130 

31 

.45 

1 

0 

0.01 

87 

11 

130 

31 

.25 

1 

0 

0.01 

73 

12 

130 

31 

0 

1 

1 

0.01 

60 

13 

130 

31 

0 

1 

0 

0.01 

60 

14 

130 

31 

.45 

0 

1 

0.025 

60 

15 

130 

31 

.45 

0 

1 

0.01 

36 

*  Calculated  speeds  are  higher  at  z  >  60m. 
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The  most  commonly  used  value  of  the  parameter        is  unity  (see,  for  example,  Ref.  3  and 

12).     However,  it  is  by  no  means  certain  that  Eqs.  13  and  19  with  "'"  ^  satisfactory 

model  of  the  radial  velocity  v^^.     Indeed,  according  to  Ref.  10  (p.  135),  v^^  is  negligible 

throughout  the  tornado  flow  except  in  the  vicinity  of  r  =  R^  (z) .     It  is,  therefore,  of 

interest  to  examine  the  influence  upon  the  calculated  missile  speeds  of  the  assumption 

V    =0.    This  can  be  done  by  assuming  k„  =  0  in  Eqs.  13  and  19. 
R  ^ 

Similarly,  while  it  is  commonly  assumed        =  1,  it  is  of  interest  to  examine  the 

case  in  which  the  value  of  the  vertical  velocity  component  v^  is  smaller  than  that  given 

by  Eqs.     17  and  23  with        =  1.     This  can  be  done,  for  example,  by  assuming  k^  =  0. 

Calculations  of  missile  trajectories  and  speeds  based  on  the  assumption  that  the 
tornado  vortex  is  described  by  Eqs.     13  through  24  were  carried  out  for  the  15  cases 
listed  in  Table  6.     In  all  cases  of  Table  6  it  was  assumed  x(0)  =  R^^,  =  46m,  y(0)  =  0, 
z(0)  =  0.     Table  6  only  includes  the  maximum  calculated  horizontal  missile  speeds  v'"^^ 
at  elevations  z  _<  60  m. 

A  discussion  will  now  be  presented  of  the  results  of  Table  6.     The  cases  corresponding 
to  the  set  of  parameters  k^  =  0.45,  k^  =  1,  k^  =  1  which,  as  previously  Indicated,  are 
most  commonly  assumed  in  missile  speed  calculations,  will  be  examined  first.     It  is 
noted  that  the  missile  speeds  for  cases  1,  5  and  9  of  Table  6  are  higher  than  the  corresponding 
speeds  of  Fig.  3,  as  shown  in  Table  7. 


Table  7  -  Comparison  between  Calculated  Speeds  Based  on  Eqs.  13  through  24 

and  on  Eqs .  5  through  9 


Table 

6 

Figure  3 

m 

Case 

max 
(m/s) 

Tnrnado 
Type 

max 
(m/s) 

0.025 

1 

90 

3* 

57 

0,025 

5 

90 

1* 

81 

0.01 

9 

91 

1* 

63 

*    See  Table  1 

This  can  be  explained  as  follows.     Assume  that        ~  '^j  ~  ^  ^'^'^  that  the  missile  starts 

from  an  initial  position  x(0)  =  R  .     As  the  missile  gains  momentum,  it  tends  to  fly  in  a 

m 

straight  line  along  the  direction  of  the  tangential  velocity,  i.e.,  at  larger  distances 

from  the  center  of  the  vortex.     In  the  case  of  a  tornado  with  a  radius  of  maximum  tangential 

speed  that  is  constant  with  height,  as  the  missile  moves  away  from  the  intial  position,  it 

will  be  subjected  to  winds  of  lower  intensity  and  its  speed  will  increase  more  slowly. 

On  the  other  hand,  if  the  radius  of  maximum  tangential  speed  increases  with  height,  and  if 

the  elevation  of  the  missile  increases  under  the  action  of  vertical  wind  speed  components 

larger  than  the  missile  weight,  then  as  the  missile  moves  away  from  the  center  and  up  it 

continues  to  fly  in  zones  of  maximum  winds  and  thus  continues  to  gain  momentum  at  a  fast 
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rate.  This  mechanism  is  modified  only  slightly  in  most,  although  not  in  all  cases,  if 
non-zero  radial  and  translation  velocities  are  present.  The  explanation  just  advanced 
is  confirmed  by  comparing  cases  4  to  1,  8  to  5,  and  12  to  9.     In  cases  4,  8  and  12,  = 

0,  i.e.,  the  radius  of  maximum  tangential  velocity  is  constant,  and  the  missile  speeds 
are  lower  than  in  cases  1,  5,  9,   in  which  that  radius  increases  with  height.     A  similar 
conclusion  is  reached  by  comparing  cases  3  to  2,  7  to  6,  11  and  13  to  10. 

The  parameter  k^,  which  controls  the  magnitude  of  the  vertical  wind  velocity  component, 
does  not  appear  to  affect  the  missile  speeds  in  a  uniform  fashion,  i.e.,  to  the  value 
=  0  there  correspond  in  certain  cases  lower  values  of  the  missile  speed  than  to  the 
value  k^  =  1  (e.g.,  case  2  vs.  case  1);  in  other  cases,  the  reverse  is  true  (e.g.,  case 
6  vs .  case  5) . 

Some  meteorologists  have  expressed  the  view  that  the  radial  velocity  component  v 

R 

is  smaller  over  most  of  the  tornado  field  than  that  described  by  either  Eq.  7  or  Eqs .  13 
and  19  with  k^  =  1[10] .     It  is  therefore  of  interest  to  examine  the  case        =  0.     It  is 
seen,  by  comparing  in  Table  6  cases  14  and  15  with  cases  5  and  9,  respectively,  that  the 
assumption  k^  =  0  results  in  considerably  lower  missile  speeds  than  the  assumption  k^  = 

1.  This  is  the  case  because,  if  k    =  0,    (i.e.,  if  v    =  0) ,  the  missile  is  ejected 

Z  K 

(i.e.,  it  flies  at  larger  distances  from  the  center  of  the  vortex  and  therefore,  in  a 

region  of  lower  speeds)  sooner  than  if        =  1,   (i.e.,  than  if  a  radial  velocity  is 

present  that  resists  the  tendency  of  the  missile  to  fly  away  from  the  center) .  This 

reasoning  is  valid  for  the  assumed  initial  condition  x(0)  =  R  ,  y(0)=  0.  it  is  conceivable 

m 

however,  that  initial  conditions  exist  which  might  result  in  higher  missile  speeds  for 
k^  =  0  than  for  k^  =  1. 

6.     AERODYNAMIC  FORCES 


In  this  work  it  was  assumed  that  the  aerodynamic  force  acting  upon  a  missile  is 

D  =  1/2  p  C^A  I  V    -  v„  I (v    -  O  (24) 
D     '     w        M  '     w  M 

where  p  =  air  density,  v^  =  wind  velocity,  v^^  =  missile  velocity  and  C^A  is  a  suitably 

chosen  quantity  with  the  dimensions  of  an  area.     As  indicated  in  Section  2,  Eq.  2  may  be 

assumed  to  be  valid  under  random  tumbling  conditions,  i.e.,  if  the  motion  is  such  that 

the  body  "presents  all  possible  aspects  to  the  flow, .    .    .  the  orientation  of  the  surface 

elements  with  respect  to  the  flow  sweeps  through  all  the  possible  angles"   (Ref .  14)  and, 

moreover,  the  angular  speed  of  the  tumbling  body  is  sufficiently  large.     In  this  connection 

it  may  be  argued  that,  instantaneously,  lift  forces  do  exist,  and  that  the  momentum  they 

impart  to  the  missile  may  not  be  negligible  if  the  rotational  motion  of  the  missile  is 

relatively  slow.     The  acceleration  of  the  missile  would  then  no  longer  be  parallel,  at 

every  instant,  to  the  vector  v    -  v    and  the  missile  trajectory  would  deviate  from  that 

w  M 

predicted  by  using  the  aerodynamic  model  implicit  in  Eq.  2.     It  is  believed,  however, 
that  the  effect  of  such  deviations  on  the  maximum  speed  of  the  tornado-borne  missile  is 
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comparable  to  the  effect  of  changes  in  the  initial  conditions  of  the  problem  such  afe  were 
stud  .led  i.n  Section  5.     For  the  purposes  of  this  report,  the  existence  of  lift  forces  which 
are  not  taken  into  account  in  the  calculations  is  believed  not  to  invalidate  the  aerodynamic 
model  used  herein. 

The  value  C^^A  that  corresponds  to  tumbling  conditions  can,  in  principle,  be  determined 
experimentally.     Unfortunately,  little  information  on  this  topic  appears  to  be  presently 
available.     Ref .  14  examines  tumbling  motions  for  supersonic  wind  conditions,  while  Ref . 
15  contains  information  on  tumbling  under  flow  conditions  corresponding  to  Mach  numbers 
0.5  to  3.5.    Hoerner  extrapolated  the  data  of  Ref.  15  to  lower  subsonic  speeds  (Ref.  16, 
p.  14-16,  Fig.  7);  according  to  this  extrapolation,  for  a  randomly  tumbling  cube  the 
quantity  C^^A  equals,  approximately,  the  average  of  the  products  of  the  projected  areas 
corresponding  to  "all  positions  statistically  possible"  by  the  respective  static  drag 
coefficients  (Ref.  16,  p.  14-16  and  P.  13-17).    An  investigation  into  this  question  is 
currently  carried  out,  within  the  framework  of  this  project,  by  Colorado  State  University 
(CSU) .    The  theoretical  and  experimental  results  of  this  investigation  will  be  reported 
in  a  separate  document  by  CSU. 

In  the  absence  of  more  experimental  information,  it  appears  reasonable  to  assume 
that  the  effective  product  Cj^A  is  given  by  the  expression 

V  =      \h  ^  S3^3>  (25) 

in  which  Cjj  A^  (i  =  1,  2,  3)  are  products  of  the  projected  areas  corresponding  to  the 

cases  in  which  the  principal  axes  of  the  body  are  parallel  to  the  vector  v    -  v.,,  by 

w  M 

the  respective  static  drag  coefficients.     In  Eq.  25,  c  is  a  coefficient  assumed  to  be 
0.50  for  planks,  rods,  pipe  and  poles  and  0.33  for  the  automobile. 
An  upper  bound  for  the  quantity  Cj^A  is  believed  to  be 

(C^A)  =  (26) 

in  which  C^^  A^  is  the  largest  of  the  quantities  C^^  A^  (i  =  1 ,  2 ,  3)  . 
The  Reynolds  number  is  defined  as  ^ 

Re  =^ 

V 

where  V  =  fluid  velocity  relative  to  the  body,  d  =  characteristic  dimension  of  the  body 

(in  the  case  of  a  cylinder,  d  =  diameter)  and  v  =  kinematic  viscosity  (v  =  1.5  x  10~^ 

m/s  for  air).    For  a  circular  cylinder  Re  -  .67  x  10^  Vd,  where  V  and  d  are  expressed  in 

m/s  and  m,  respectively.    For  Re  >  4  x  10^,  i.e.,  the  Reynolds  number  is  in  the  supercritical 

range  and  it  may  therefore  be  assumed,  conservatively,  Cj^    =  0.7  (see  Ref.  8,  p.  67).  In 

the  case  of  the  1  inch  (2.54  cm)  rod,  it  may  be  assumed  tftat  Re  is  in  the  subcritical 

range  even  for  velocities  V  of  the  order  of  100  m/s  and,  therefore,  that  C      =    1.2  (Ref.  8). 

°1 
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7.     SPEEDS  OF  SELECTED  POTENTIAL  TORNADO-BORNE  MISSILES 


In  this  section  calculated  speeds  of  selected  potential  tornado-borne  missiles  will 
be  given,  based  on  the  following  assumptions: 

(1)  The  model  of  the  tornado  vortex  consisting  of  Eqs.  5  through  9  is  valid,  with 
the  parameter  values  corresponding  to  cases  1,  2  and  3  of  Table  1. 

(2)  The  initial  conditions  are  x(0)  =  R  ,  y(0)  =  0,  z(0)  =  40m  (for  comments  on 

m 

the  initial  condition  z(0)  =  40m,  see  p.  lo  of  this  report);  v^(0)  =  '^My^^^  "  ^Mz^^^ 
=  0.    Assumptions  (1)  and  (2)  just  described  were  used  in  calculating  the  values  on 
curves  1,  2  and  3  of  Fig.  3. 

(3)  The  effective  product  Cj^A  is  given  by  Eq.  25. 

The  results  of  the  calculations  based  on  these  assumptions  are  shown  in  Table  8. 

The  missile  speeds  of  Table  8  are  based  on  a  set  of  assumptions  which,  while  reasonable, 
might  in  certain  cases  not  correctly  reflect  the  actual  physical  phenomenon.     It  follows 
from  Sections  5  and  6  that  the  order  of  magnitude  of  uncertainties  in  the  estimates 
of  maximum  missile  speeds  can  in  certain  cases  be  as  high  as  50%  or  even  more.  Whether 
or  not  actual  missile  speeds  will  be  higher  than  those  listed  in  Table  8  depends  in 
large  measure  on  the  extent  to  which  the  tornado  flow  model  consisting  of  Eqs.  5  through 
9  (the  so-called  Rankine  vortex)  is  realistic.     In  particular,  if,  as  suggested  in  Ref . 
10,  the  radial  and  vertical  velocity  components  in  a  tornado  are  actually  lower  than 
those  given  by  this  model,  it  could  be  expected — all  other  conditions  being  equal — 
that  the  predictions  of  Table  8  are  conservative.     If,  on  the  other  hand,  the  actual 
tornado  flow  is  more  closely  represented  by  Eqs.  13  through  24  with  certain  unfavorable 
values  of  the  parameters  included  in  these  equations,  then  higher  missile  speeds  than 
those  of  Table  8  may  occur,  as  shown  in  Tables  6  and  7. 

The  speeds  of  Table  8  may  also  be  exceeded  if  unfavorable  initial  conditions  obtain. 
The  uncertainties  pertaining  to  the  tornado  flow  modeling  are  due  to  the  lack  of  reliable 
information;  on  the  other  hand,  those  pertaining  to  the  initial  conditions  are  a  consequence 
of  the  probabilistic  nature  of  the  problem.     Probabilities  of  occurrence  may  be  assigned 
to  each  set  of  initial  conditions.     The  probability  that  (a)  the  wind  speed  will  reach 
the  intensity  levels  of  Table  1,   (b)  that  a  missile  starts  from  a  highly  unfavorable 
set  of  initial  conditions  and  (c)  that  it  will  hit  a  certain  installation  with  a  speed  V 

H 

can  be. expected  to  be  negligibly  low.     Such  probabilities  can  in  principle  be  evaluated 
using,  for  example,  procedures  similar  to  those  outlined  in  Ref.  18.    Attempts  to  calculate 
such  probabilities  are  beyond  the  scope  of  this  work. 
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The  uncertainties  regarding  the  actual  aerodynamic  drag  coefficients  constitute 

another  source  of  error.     It  is  noted  that  the  curves  in  Fig.  3  are  S-shaped.  Large 

errors  in  the  assumed  value  of  the  quantity  C^^A  may  therefore  in  certain  cases  result  in 

considerable  changes  of  the  estimated  value  of  V™t^    For  example,  if,  for  a  12  in  Sch. 

2 

40  pipe  (entry  6  in  Table  8)  it  was  assumed  C^A/m  =  0.0066    kg/m    instead  of  C^A/m  =  0.0033 

kg/m^,  it  follows  from  Fig.  3  that  V^^^  =  40  m/s,  rather  than  V™^^  =  7  m/s,  as  in  Table 

8.     It  is  interesting  to  note,  on  the  other  hand,  that  as  long  as  a  change  in  the  assumed 

value  of  Cj^A/m  does  not  displace  the  point  from  the  upper  or  from  the  lower  branch  of  an 

S-shaped  curve,  the  sensitivity  of  V^^^    to  even  considerable  changes  in  the  value  of    C  A/m 

is  fairly  small.     For  example,  if  for  the  6  in  Sch.  40  pipe  (entry  2  in  Table  8)  it  is  assumed 

that  C  A/m  =  0.0086,  then,  for  a  tornado  type  1,  V™^^  =  61  m/s  (see  Fig.  3),  whereas  to 

ax 

the  assumption  C^A/m  =  0.0043  there  corresponds  V?  ^=  52  m/s,  i.e.,  to  an  error  of  100%  in 
D  H 

the  value  of  C  A/m  there  corresponds  an  error  of  only  17%  in  the  estimated  value  of  V 
D  n 

Finally,  it  must  be  noted  that  the  actual  missiles  may  have  properties  that  are 
more  unfavorable  than  those  listed  in  Table  8;  for  example,  the  case  is  mentioned  in 
Ref .  19  of  a  beam  attached,  during  its  flight,  to  a  portion  of  a  carport  roof  which 
considerably  increased  the  surface  area  of  the  missile  and,  therefore,  the  parameter 
CjjA/m. 

8.  CONCLUSIONS 

In  the  preceding  section  calculated  maximum  speeds  of  tornado-borne  missiles  are 
given  in  Table  8,  based  on  a  set  of  assumptions  believed  to  be  reasonable.     However,  in 
assessing  these  speeds,  it  must  be  recognized  that: 

1.  The  problem  of  determining  tornado-borne  missile  speeds  has  a  probabilistic 
character.    As  shown  in  the  body  of  the  report,  unfavorable  initial  conditions  may 
obtain — to  which  there  correspond  relatively  low  probabilities  of  occurrence — for 
which  the  maximum  missile  speeds  would  be  higher  than  in  Table  8.     The  estimation  of 
such  probabilities  is  beyond  the  scope  of  the  investigation  covered  by  this  report. 

2.  Estimates  of  tornado-borne  missile  speeds  are  also  affected  by  significant 
uncertainties  with  regard  to:     (a)  the  detailed  structure  of  the  tornado  flow  and 
(b)  to  the  aeordynamic  behavior  of  the  missile.     Under  certain  assumptions  regarding 
one  or  both  of  these  factors,  calculated  missile  speeds  can  be  higher  than  those 

of  Table  8.     However,  it  is  believed  that  the  assumption  used  to  derive  the  values 

of  Table  8  are  conservative.     In  particular,  it  is  believed  that  the  actual  vertical 
wind  speeds  are  lower  than  indicated  by  Eq.  9,  so  that  the  relatively  heavy  missiles 
would  tend  to  hit  the  ground  sooner  than  calculated  on  the  basis  of  this  equation, 
with  a  consequent  reduction  in  the  calculated  maximum  missile  speed. 

In  spite  of  the  many  uncertainties  involved,  the  writers  believe  that  the  assumptions 
used  to  estimate  the  speeds  of  Table  8  are  sufficiently  conservative  for  purposes  of 
nuclear  power  plant  design.     It  is  the  writers'  judgement  that,  although  higher  values  of 
tornado-borne  missile  speeds  are  conceivable,  their  probabilities  of  occurrence,  for  any 
given  tornado  strike,  are  low. 
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APPENDIX  A  -  EQUATIONS  OF  MOTION 

Al.l    Absolute  and  Translating  Frames  of  Reference 

Consider  the  absolute  frame  of  reference  Oj^  xyz  (Fig.  Al.l)  and  a  frame  0^  x'y'z 

such  that  0^2  is  parallel  to  Oj^z.     It  is  assumed  that  the  frame  O2  x'y'z  translates  with 

respect  to  0^^  xyz  with  a  constant  velocity,  the  components  of  which  (in  the  Oj^  xyz  frame) 

are  (v„  ,        ,  0) .    The  angle  between  0-x'  and  O.x    being  denoted  by  6  , 
Tx      Tv  2  1  c 


""tx 

®c  -  ,2^2  ,1/2  (Al.l) 
(^x  ^  ^y> 


sin  0 


Ty  

,2.2  .1/2 


(A1.2) 


A1.2    Vectors  of  Position 

Using  the  notations  of  Fig.  Al.l,  the  vectors  of  position  of  the  particle  P  in  the 
absolute  and  in  the  translating  frame  of  reference  may  be  written  as  follows: 

Absolute  frame  (with  origin  0^^) 

R=xT+yT   +    zk    =    r  cos  0  T   +  r  sin  0  j    +    z  k  (A1.3) 
where  T,  J,  k  are  unit  vectors  along  the  axes  O^x,  O^^y,  O^^z  respectively. 
Translating  frame  (with  origin  O2) 

R'    =    x'  i'    +    y'  J'    +    z  k    -    r'  cos  0'  i'    +    r'  sin  0'  J'  +  z  k  (A1.4) 

where  1?,  j',  k  are  unit  vectors  along  the  axes  O2X' ,  02y',  02Z»  respectively 

A1.3    Transformations  of  Unit  Vectors 

The  unit  vectors  r,  0,  k  of  the  revolving  frame  of  reference  (see  Fig.  Al.l)  may 
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be  written  in  terms  of  the  unit  vectors  of  the  translating  frame  02X'y'z  and  the  angle 
between  the  lines  O^P  and  O^x",  where  O^x"  is  parallel  to  O2X' ,  as  follows: 


r    =    cos  e'  i'  +    sin  e'  j' 


The 


e  =  -sin  6'  i'  +  cos  e'  j'  (A1.5) 
k   =  k 

imit  vectors  i' ,  j',  k,  of  the  translating  frame  0^  x'y'z    may  be  written  in 


terms  of  the  unit  vectors  of  the  absolute  frame  of  reference  and  of  the  angle  6^  as 
follows: 


i     =    cos  e      i    +    sin  6  j 
c  c*^ 


j'    -  -sin  e    i    +    cos  e    j  (A1.6) 
c  c 


A1.4    Expression  of  Wind  Velocity  in  an  Absolute  Frame  of  Reference 

Let  the  vortex  wind  velocity  with  respect  to  the  translating  frame  of  reference  be 
denoted  v^,  and  let  its  components  along  the  vectors  r,  9  and  k  be  denoted  by  v^^.  v^g 

and  v^j^,  respectively,  i.e.. 


V     "    V       r    +    V  „  9     +    v  ,  k  (A1.7a) 
(0  ur  u)6  (i)k 


With  respect  to  the  absolute  frame  of  reference,  v^  may  be  written  as 


v-vl    +    vj+vk  (A1.7b) 
u  ux  coy  u)Z 


The  components  v    ,  v     and  v     are  obtained  by  substituting  Eqs.  A1.5  and  A1.6  into 
ux     uy  (DZ 


Eq.  A1.7a,  i.e.,  in  matrix  form. 
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ux 


0)  z 


cos    e    -  sin  e  0 
c  c 


sin    e       cos  e  0 
c  c 


cos    e'     -sin    e'  0 


lin    e'      COS    e'  0 


V 

V 

0)6 

(A1.8) 


The  total  wind  velocity  consists  of  the  sum  of  the  vortek  wind  velocity  and  the 
translating  velocity  of  the  vortex,  i.e., 


V  i  +  v  J+v,  k 
wx  wy  wk 


(A1.9) 


where 


V  ■  V  +  v_ 
wx       ojx  Tx 


V      "V      +  v_ 
wy         uy  Ty 


(Al.lO) 


V        =  V 
WZ  U)Z 


and  V    ,  V    ,  v     are  given  by  Eqs.  A1.8.    In  Eqs.  Al.lO  and  A1.8,  the  quantities 

v_  ,  v_  ,  V    .  V       V  ,   and  6  are  specified.    The  quantity  9  '  is  a  function  of  time  and 
ix      iy      lor      0)0      a)K  c 

is  obtained  from  the  relations 


cos  0      =  — 


sin  e'  =  ^ 
r' 


(Al.ll) 


,  ,  ,2  ^  ,2.1/2 y 

r '  «  (x'    +  y'  ) 
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in  which  x'  and  y'  are  determined  as  explained  in  the  following. 
Let  the  initial  conditions  of  the  problem  be 


t    =    initial  time 
o 

(x  ,  y  ,  z  )      =      coordinates  of  particle  at  time  t    in  the  absolute  frame 
o     o     o  ,  o 

(v      .         ,  v„    )  =    velocity  components  of  particle  at  time  t    in  the 

^.PXO'      Pyo'      PZO  J  r  r  q 

absolute  frame 

(x02q»  y^2o*        ~  P°^^^^°°  (oJ^^Sin  of  the  translating  frame)  at  time 

t  ,  in  the  absolute  frame, 
o 

At  time  t,  the  position  vector  of  the  origin  O2  is  ?  , 


=  [xO^^  +        (t  -  t^)]  i  +  [y02^  +        (t  -t^)  j] 


BO  that  for  any  time  t  >_  t^  we  have 


(A1.12) 


R'  =  R  -  RO, 


(A1.13) 


where 


R=xi+yJ  +  zk,  i.e.. 


^R' 


=R' 


xO_      +    v„    (t  -  t  ) 

2o    ,       Tx  p 


yO„  +  v_.  (t  -  t  ) 
■'20       Ty  o 


(Al.lA) 


Since 


Xr.  i    +    Yr.  j    +    Zr.  k    =  x'  i'    +    y'  j'    +  z'  k 


(A1.15) 


it  follows,  using  the  inverse  of  the  transformation  A1.6, • 


L  z'  J        L  0 


cos  6      sin  6  0 
c  c 


-sin  e     cos  e  0 
c  c 


0 


^R- 


«-  ^R' 
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(A1.16) 


where  x^^,,  y^,,  Zj^i,  are  given  by  Eq.  Al.14.    The  quantities  x',  y',  z',  in  Eq.  Al.ll,  and 
thus  the  quantity  6*  in  Eqa  Al.lO,  A1.8  are  expressed  in  terms  of  the  coordinates  of  point 
f  with  respect  to  the  absolute  frame,  and  of  the  quantities  xO^^,  y02o*  ^xx*  ^Ty*  ^o  ^' 
The  problem  of  expressing  the  components  In  the  absolute  frame  of  the  wind  velocity  at  point 
P  (x,  y,  z)  in  terms  of  the  vortex  wind  velocity  components  v^^,  v^^,  v^j^,  of  the  translation 
velocity  components  v^^,  v^^,  of  the  initial  conditions  and  of  the  coordinates  x,  y,  z  of 
the  point  P  Is  thus  solved. 


A1.5    Equations  of  Motion 

The  equations  of  motion  nay  be  trriccen  as 


1/2  P  V!v^Jv^j-gk 


(A1.17) 


where  m  •■  mass  of  particle,  p  ■  air  density,  C^^  "  drag  coefficient,  A  >  area  of  the  particle, 

•  •  • 

Vj^gj^  -  (v^^  -  x)  1   +  (v     -  y)  T  •♦•  <v^^  -  z)  k,  and  g  -  acceleration  of  gravity. 


APPENDIX  B  -  ANALYTIC  SOLUTION  TO  THE  EQUATIONS  OF  MOTION  FOR  A  UNIFORM  WIND  FIELD. 

In  the  case  of  a  uniform  wind  field  an  analytic  solution  of  the  equations  of  motion 
of  a  particle  may  be  obtained  as  shown  herein.    This  solution  was  used  to  test  the  computer 
program, in  which  the  wind  field  subroutine  was  suitably  modified.    The  following  assumptions 
were  used : 

1)  The  initial  velocity  of  the  particle  Is  zero 

2)  The  notion  occurs  In  the  x  -  z  plane  only 

3)  The  wind  velocity  vector  Is  at  all  points  parallel  to  the  horizontal  axis  0^  and 
has  the  constant  magnitude  v  . 
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The  equation  of  motion  in  the  horizontal  direction  is 


where 


.2  ,  „ 

dx  ,         dx^  2  .  / . o  T  \ 

—1  =  "  ^^w  -  (^2.1) 
dt 


^  =  I  CjjA,  or 


f     '^(f)         =     /•  a  dt  + 


Cj^  '  (A2.2) 


With  the  change  of  variables  u  =       -  v^,  du  =  d         Eq  A22  becomes 


J  ^   =  a  t  +  (A2.3) 


-    =  a  t  +  (A2.4) 

dt  " 

Integrating  Eq.  A2.5,  there  follows 

X  =  V  t  -  /   +  (A2.6) 

w     y  at  +  2 

With  the  change  of  variables  a  t  +       =  u,  a  dt  =  du,  dt  =         Eq.  A2.6  becomes 

x  =  vt--/—   +  (A2.7) 

w       a  /  -  u  2 

X  =  V  t  -  -   In  (at  +  C.)  +  C„  (A2.8) 
w       a  1  2 

The  constants  of  integration  are  determined  as  follows : 

At  t  =  o ,  X  =  X     and  4^    =  o ,  i . e . , 

o  dt 


27 


-  -  -  In  C-  +  C- 
o         a         1  2 


(from  Eq.  A2.8) 


(A2.9) 


1 

^w- c: 


(from  Eq.  A2.5) 


(A2.10) 


Thus 


-    7"  (A2.11) 

V 


C-  -  X  In  V 

2       o      a  w 


(A2.12) 


and 


x"x    +  vt  ln(avt+l) 

o       w       a  w 


(A2.X3) 


28 


APPENDIX  C 


DOCUMENTATION, 


AND 


SAMPLE  INPUT  AND  OUTPUT 


OF  COMPUTER  PROGRAM 


29 


TnHNArO**  IMM  1  )  .*«A  |N«CX  (  1  ) 


1 

C 

2 

c 

NATIONAL   BUREAU   OF  STANOAKDS 

3 

c 

APPLIFO    MATHF.MATICS  DIVISION 

4 

c 

5 

c 

TOWWDl 

6 

c 

7 

c 

A   FCRTHAN   APPLICATIONS    PROGRAM   TO   STUOV    THF   MUTfGN  OF   AN  OBJFCT 

B 

c 

UNOCP    THE    INFLOFNCr    OF    A    TOPNADO    W I  NO  FIFLO, 

q 

c 

1  0 

c 

MAWTIN  COROES 

1 1 

c 

1  2 

c 

APRIL  197b 

I  3 

c 

1  4 

1  5 

c 

16 

c 

PU«POSF. 

1  f 

c 

I  a 

c 

THIS  FCRTPAN   PROGRAM    INTFGRATFS    IN    T I MF    THF    FOLLOWING   EQUATIONS  OF 

1  Q 

c 

MCTICN   FOR    A   RIGID    f^OOy    ACTED  ON    RY    A   DRAG    FORCE    IN  A 

20 

c 

GRAVITATIONAL  FIFLC* 

2  1 

c 

22 

c 

0**2(PX»   /    01**?   =    ALPHA   *  VRX 

23 

c 

24 

c 

0**2(py)    /    DT**2    =    ALPHA    *  VRY 

25 

c 

• 

26 

c. 

0**2(PZ»    /    OT*«2  =    ALPHA    *    VRZ    -  G 

27 

c 

28 

c 

ANO 

29 

c 

30 

c 

Alpha  =  cd  •  a  *  rho  *  vr  /  2  • 

31 

c. 

32 

c 

33 

c 

34 

c 

0**2(X)    /    0T»*2         IS    THF    SECOND  DERIVATIVE    OP    X    WITH  RESPECT 

35 

c 

TO  TIME. 

36 

c 

37 

c 

PX                                       ARE    THE    X.    Yt    2    COORDINATES   OF  THE 

38 

c 

PV                                         PARTICLE    AS    A   FUNCTION   OF    T 1 MF    IN  THE 

39 

c 

P2                                       ABSOLUTE  FRAME, 

40 

c 

4  1 

c 

VRX                                      ARE    THE    X.    Y,    Z    CUMFCNtNTS    OF    THE  RELATIVE 

42 

c 

VRY                                     VELOCITY   OF    THF    WIND    WITH  RESPECT    TO  A 

43 

c 

VRZ                                     STATIONARY    PARTICLE    IN   THF    ABSOLUTE  FRAME. 

44 

c 

45 

c 

G                                            IS    THE    MAGNITUDE    OK  GRAVITY. 

46 

c 

4"' 

c 

CO                                         IS    THE    DRAG  COEFFICIENT. 

4a 

c 

4  9 

c 

A                                          IS    THE    AREA    ASSOCIATED    *ITH   THE  BODY. 

60 

c 

51 

c 

RHO                                       IS    THE    DENSITY    OF  AIR. 

62 

c 

53 

c 

VR                                       IS    SORT (VRX    **    2    ♦    VRY    »♦   2    ♦    VRZ   **  2). 

54 

c 

55 

c 

THIS    SYSTEM   CF   THREE.    SECOND    ORDER    EQUATIONS    IS    CONVERTED   TO  A 

56 

c 

SYSTEM    OF    SIX.    FIRST    ORDER    EQUATIONS   BY    A    STANDARD  TECHNIQUE. 

5-» 

c 

THE   rOUfVALENT    FIRST    ORDER   SYSTEM    IS    INTEGRATED   RY    THE  ODC 

30 


58 

C 

SOLVP«.                                            .  . 

5S 

c 

60 

c 

THE   COMPONPNTS   OF    THE    PARTICLE    VELOCITY    AND    THE    WIND   VELnCITY  APE 

6  1 

c 

*>ITH   RESPECT    TO   AN   ABSOLUTE   FRAME,    ThF    WIND    VELOCITY    IN  THE 

e2 

c 

ABSOLUTE   FRAME    IS    THE   SUM   OF    THE    WIND   VELOCITY    FRQM   A  STATIONARY 

6  3 

c 

TORNADO    VORTEX    AND    THE    TRANSLATION    VELOCITY    OF    THE    TORNAOC  VGfiTFX. 

64 

c 

65 

c 

INITIAL  CONDITIONS, 

66 

c 

, 

67 

c 

AT    TIME   =   TO    THE   X,    Y,    Z    COORDINATES    AND   VELOCITY  COMPONENTS 

6H 

c 

OF    THE  PARTICLE    ARE    SPECIFIED,    IN  ADOITICN,    THE   X,  Y 

69 

c 

COORDINATES    OF    THE   ORIGIN   OF    TmE    TRANSLATING  FRAME    WHICH  THE 

70 

c 

TORNADO    IS    STATIONARY    IN    ARE  SPECIFIED, 

71 

c 

72 

c 

CONDITIONS    INDEPENDENT   OF  TIME, 

77 

c 

74 

c 

THE    X,    Y    VELOCITY    COMPONENTS    OF   THE    TRANSLATING  FRAME  ARE 

75 

c 

CONSTANT. 

76 

c 

77 

c 

FOR  FURTHER    DETAILS  CONSULT    REFERENCE  2), 

78 

c 

79 

c 

APPLICABILITY    ANC   RFSTr  ICT  ICNS . 

80 

c 

81 

c 

THIS    PROGRAM    IS    INTENDED   FOR   WORK    THAT  REQUIRES   MODERATE  ACCURACY 

P2 

c 

AND    MODERATE   TIME    PERIODS    OF  INTEGRATION, 

83 

c 

8A 

c 

THE   CURRENT    VERSION   PRODUCES   PLOTS  OF   AT   MOST    101    DATA    POINTS.  THE 

85 

c 

DATA    IS  FROM    THE    TABULAR  RESULTS   GENERATED    AT   THE   FIRST    10  1  PRINT 

86 

c 

STEPS.    TO  PREVENT   A    TRUNCATED  PLOT   FROM   BEING  GENERATED    RE  SURE 

37 

c 

THAT    TIMEIT    AND   FXTMIT    SATISFY    T  IME  IT    ,LT ,    100    ♦  FXTMIT,  HLPE 

88 

c 

TIMEIT    IS    THE    TIME    INTERVAL    OF    INTEGRATION    AND   FXTMIT    IS  THE  TIME 

89 

c 

INTERVAL    BETWEEN  SUCCESSIVE    PRINT  STEPS. 

90 

c 

91 

c 

PROGRAM  STRUCTURE. 

92 

c 

93 

c 

SUBPROGRAM  DIRECTORY. 

94 

c 

95 

c 

COMMENT      ADDITIONAL    SUBPROGRAMS   OTHER   THAN  FORTRAN  INSTRINSIC 

96 

c 

FUNCTIONS    AND  FORTRAN   BASIC   EXTERNAL  FUNCTIONS 

97 

c 

NECESSARY   FOR    THE   EXECUTION   CF    THIS    MAIN  PROGRAM  ARE 

98 

c 

LISTED   BELOW    IN    ALPHABETICAL    ORDER.    EACH  IS 

99 

c 

ACCOMPANIED   BY    A   SHORT    FUNCTIONAL  DESCRIPTION. 

100 

c 

101 

c 

OFUN             THIS   SUBROUTINE  COMPUTES   THE  RIGHT   HAND    SIDE   OF  THE 

102 

c 

FIRST    ORDER    CDE    SYSTEM.    DY    IS   THE  LOCAL    NAME  OF  OFUN 

103 

c 

IN  VOAOAM. 

1  04 

c 

105 

c 

DRAG             (USER    SUPPLIED.)    THIS   SUBROUTINE  COMPUTES   THE  DRAG 

106 

c 

COEFFICIENT. 

107 

c 

108 

c 

INPUT           THIS    SUBROUTINE    READS    IN   OR   CONTROLS   THF    READING    IN  OF 

109 

c 

ALL  RELEVANT  DATA    AND  PARAMETERS  EXCEPT    THE  PARAMETERS 

110 

c 

REQUIRED   BY    THE   ODE   SOLVER   AND    WHICH   AFFECT  OUTPUT, 

111 

c 

112 

c 

OUTPUT        THIS   SUBROUTINE    READS    IN  PARAMETERS   THAT   CONTROL  WHICH 

113 

c 

PLOTS   ARE    TO   BE    PRODUCED   AND  SUPERVISES   THE  PRODUCTION 

1  14 

c 

AF   ALL  OUTPUT, 

115 

c 

31 


I  16 

C 

I  I 

c 

1 1  e 

c 

I  19 

c 

120 

c 

12  1 

c 

122 

c 

123 

c 

124 

c 

125 

c 

126 

c 

127 

c 

12« 

c 

129 

c 

1  30 

c 

1  3  1 

c 

1  32 

c 

133 

c 

134 

c 

135 

c 

1  36 

c 

1  37 

c 

1  3fl 

c 

1  39 

c 

140 

c 

1  4  1 

c 

142 

c 

143 

c 

14  4 

c 

145 

c 

146 

c 

147 

c 

148 

c 

149 

c 

1  50 

c 

151 

c 

152 

c 

153 

c 

154 

c 

155 

c 

156 

c 

157 

c 

15a 

c 

159 

c 

160 

c 

161 

c 

162 

c 

163 

c 

164 

c 

165 

c 

166 

c 

167 

c 

168 

c 

169 

c 

170 

c 

171 

c 

1  72 

c 

17  3 

c 

PLOT  THlf^   SUHROUTINE    PRODUCES    A   (JNt    PAGE   PRINTP.P  PLOT. 

TCnWOF         (USER   SUPPLIED.!    THIS   SUBHOUTINF   COMPUTES   THE  PAOIAL. 

ANGULAR.  AND  Z  COMPONrNTS  DF  THE  *IND  AT  THE  PAWT ICLF 
AS   MEASURED    IN    THE    TRANSLATING  FRAME. 

THIGFS         THIS    SLDMUUTINF    COMPUTES    THE    COSINE    AND    SINE   OF  THE 

BASE  ANGLE  OF  A  RIGHT  TWIANGIE  ftlTH  BASE  AND  VERTICAL 
SIDE  GIVEN. 

VCADAM        THIS   SUHHOUT  INE    IS    T  HF    ODF    SOLVtR.    EACH   CALL   DCES  ONE 
INTEGRATION  STEP. 

EXeCLTION  TREE. 


COMMENT 


THE   EXECUTION   TRFF    IS   A    STRUCTURE    THAT  OFSCRIRES  THE 
SUBPROGRAM    INTERACTION    DURING   EXECUTION.    THE    TREE  IS 
COMPOSED  OF    A  NUMBER    OF    LEVELS.    EACH   LEVEL    IS  CCMPOSED 
OF   ONE    OR  MOPE   BL'JCKS.    A   BLOCK    IS  COMPOSED   OF  A 
CONTIGUOUS    SET    OF    TREE    ELEMENTS.    FACH   TRFF  ELEMENT 
HAS   THE  FORM 

NUMHIINAMF.    KUMS2 ) 

WHERE 

NUMB  I  IS   THE    NUMBER    OF    THE   ELEMENT    IN  THE 

CURRENT  LEVEL. 

NAME  IS    THE    NAME    OF    ThF    SUBPROGRAM  ASSOCIATED 

WITH   THAT    ELEMENT , 

NUMB2  IS   THE   NUMBER  OF    THE   FIRST   ELEMENT  OF 

THE    BLOCK    AT    THE    NEXT    GREATER  LFVEL  WHICH 
CONTAINS   ALL    SU^P'^OGPAMS  CALLED  DV 
SUBPROGRAM    NAME.    IF   NUMH2    IS    ZERO  THFN 
SUBPROGRAM    NAME    CALLS    NO  SUBPROGRAMS. 


LEVEL  0 
KMAIN      •  I) 


LEVFL  1 

ic INPUT  •  n 

2(VOADAM,  4) 
3  (OUTPUT.  5) 


LEVEL  2 

KTORWDF,  0» 
2 (DRAG  .  0) 
3(TRIGFS.  0) 


LEVEL  3 

U TORWDF .  0) 
2  (  OR  A  G  .  0  ) 
3(TRIGFS,  0) 


4(0FUN  . 

5(  TOR  WDF  , 
6(CRAG  . 
7  {  TRI  GFS, 
8(PL0T  , 


1  ) 

0) 
0) 
0  » 
0  > 


PROGRAM  I/O. 
I NPUT. 

ROUT INE. 
MA  IN 
INPUT 


LOGICAL  UN  IT  (S  )  . 

5 
5 


EXECUTION  TREE  ELEMENT 
LFVEL  0,  ELFMLNT  1 
LE  VEL      1  «        ELEMENT  I 
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174 

C 

l-'S 

C 

OUTPUT 

5                                  LFVFL  1. 

ELEMENT  3 

1  76 

C 

17-r 

c 

TOR*DF 

5                                  LEVLL  2, 

ELEMENT  I 

178 

c 

17<5 

c 

DRAG 

5                                  LEVEL  2, 

ELEMENT  2 

leo 

c 

1  ei 

c 

OUTPUT, 

182 

c 

183 

c 

ROUTINE. 

LOGICAL  UNIT(S),  EXECUTION 

TREE  ELEMENT. 

164 

c 

185 

c 

HA  IN 

6                                  LEVEL  0, 

ELEMENT  1 

1  86 

c 

187 

c 

OUTPUT 

6                                  LE  VEL      I  , 

E  LE  ME  NT  3 

lea 

c 

189 

c 

TCBWDF 

6                                LEVEL  2. 

ELEMENT  5 

190 

c 

1<31 

c 

DRAG 

6                                  LEVEL  2. 

ELEMENT  6 

192 

c 

193 

c 

PLOT 

6                                  LEVFL  2. 

ELEMENT  8 

194 

c 

1  95 

c 

PROGRAM    USAGE  (INPUT). 

196 

c 

197 

c 

DATA  GROUPINGS. 

198 

c 

199 

c 

GROUP   1 . 

200 

c 

20  1 

c 

SUBGROUP 

1.1.      TORNADO   WIND  FIELD  PARAMETER 

CARD( S) . 

202 

c 

203 

c 

CARD. 

VARIABLE  LIST. 

F  ORMAT  . 

204 

c 

205 

c 

SEE 

SUBROUTINE    TORWDF    FCR   ThE   EXACT    SET    OF  PARAMETERS 

206 

c 

TO 

READ    IN   AND    THEIR  FORMAT. 

207 

c 

208 

c 

SUBGROUP 

1.2.      DRAG   COEFFICIENT   PARAMETER  CARD(S). 

209 

c 

210 

c 

SEE 

SUBROUTINE   DRAG  FOR    THE   EXACT  SET 

OF    PARAMETERS  TO 

211 

c 

RE  AC    IK   AND   The  IP  FORMAT. 

212 

c 

213 

c 

GROUP  2. 

2  14 

c 

215 

c 

SUBGROUP 

2.1.      PARTICLE   PARAMETER  CARD. 

216 

c 

217 

c 

CARD. 

VARIABLE  LIST. 

FORM  AT . 

216 

c 

2  19 

c 

1 

AREA.  MASS 

2C  12  .0 

220 

c 

221 

c 

VARIABLE 

DEFINITION   TAHLF    (GROUP  2). 

222 

c 

223 

c 

AREA 

ARE   THE    AREA.    MASS   UF    THE  PARTICLE. 

224 

c 

MASS 

225 

c 

226 

c 

GROUP  3. 

227 

c 

228 

c 

SUBGROUP 

3.1.       INITIAL   CONDITION  CARDS. 

229 

c 

230 

c 

CARD. 

VARIABLE  LIST. 

FORMAT . 

231 

c 
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23?  C  I  XPO,     YPO ,    ZPO,    VXPO«    VYPO,    VZPO  hF\?,0 

233  C 

234  C  2  XOTFt    YOfF  ,    VXTF,    VYTF,    TO  5F12.0 

?3S  C 

23ft  C  SU8GH0UP    3.?.      FINAL    CONDITION  CARD. 

2  3T  C 

23R  C  3  TIMFIT  E12.0 

23q  c 

240  C  SUBGROUP    3,1.      OOE    SULVFP    PARAMETER  CAPO. 

241  C   '     :  . 

242  C  .  4  HI.    FXTMIT,    EPSI  3ei2.0 

243  C 

244  C  SUBGROUP   3.4.      PLOT    AND    VALIDATION   PrtAMETTP  CARO. 

245  C 

246  C  .5  WTPT2Y,    WTPTYX,    WTPT7X,    WTPTXT.  712 
J47  C  WTPTSX.    WTPTHT.  WTVALO 

248  C 

249  C  VARIABLE   DEFTMTICN   TABLE    (GROUP  3). 

250  C 

251  C  .        ,  XPO  ARE   T»-E    INITIAL    X.    Y,    Z  CC.OROINATES    HE  THE 

252  C  YPO  PARTICLE. 

253  C  ZPO 

254  C 

255  C  VXPO  APE    THE    INITIAL    X»    Y,    Z   VELOCITY    COMPONENTS  OF 

256  C  VYPO  THE  PARTICLE. 

257  C  ■  VZPO 

258  C 

259  C  XOTF  ARE    THE    INITIAL    X,    Y    COORDINATES    OF    THE  ORIGIN 

260  C  YOTF  OF    THE    TRANSLATING   FRAME    *HICH    THE    TORNADO  IS 

261  C  STATICNARY  IN. 
26?  C 

263  C  VXTF  ARE   THE   CONSTANT    X.    V   VELOCITY  COMPONENTS   OF  THF 

264  C  VYTF  TRANSLATING  FRAME. 

26  5  C 

266  C  .  TO  IS    THE    INITIAL  TIME. 

26-^  C 

268  C  TIMEIT         IS   THE    INTERVAL    OF  INTEGRATION. 

269  C 

270  C  HI  IS    THF    INITIAL    TIME    STEP,    A   GOOD    VALUE    TO  USE 

271  C  IS    I .OF-4 

272  C 

273  C  FXTMIT         CONSTRAINS    THE    ODE    SOLVER    TO    PRODUCE    RESULTS  FQi^ 

274  C  PRINTING   AT    EQUALLY    SPACED   TIME    INTERVALS  OF 

275  C    ■  THIS  SIZE. 

276  C 

277  C  EPSI  IS   THE   LOCAL    ERROR   TOLERANCE  USFO   BY    THE  ODE 

278  C  SOLVER.    A    GOOD    VALUE    TO    USE    IS  l.OE-4. 

279  C 

280  C  THE  NEXT  6  VARIABLES  ARE  FLAGS  FOR  PRODUCING  PLOTS.  IF 
28  1  C  A   FLAG    IS    SET    TO    1     THE    CORRESPONDING   PLOT    IS  PRODUCED. 

282  C  OTHERWISE.    THE   PLOT    IS   NOT  PRODUCED. 

283  C 

284  C  VKTPTZV         IS    THE   FLAG    FOR    A    PLOT    CF    2   VERSES  V. 

285  C 

286  C  VkTPTYX         IS    THE    FLAG   FOR    A    PLOT    OF    Y   VERSES  X. 

287  C 

288  C  WTPTZX         IS   THE  FLAG  FOR    A   PLOT   OF    Z    VERSES  X. 

289  C 


34 


29  0 
2Q1 
29  2 
293 

29  4 
295 
2Se 
297 
298 
299 
300 

30  1 
202 
303 
304 
3  05 

2oe 

30-^ 
308 
309 

31  0 
31  I 
3  12 
313 
314 
315 
316 
317 
318 
319 
320 
321 
322 
323 
324 
325 
226 
327 
328 
229 
330 
231 
332 
3?  3 
334 
335 
236 
337 
338 
339 
34.0 
341 
342 
34  3 
344 
245 
346 
34  7 


C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c. 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


WTPTXT  IS    THE    FLAG   FOR  A    PLOT    OF    X    Vfc'RSFS  TIME. 

WTPTSX         IS    THE   FLAG  FDR  A    PLOT    QF    SPEED    VERSES  X. 

(KTPTHT         IS    THE   FLAG   FOR  A    FLCT    OF    HCRIZONTAL  SPFED 
VFRSFS  TIME. 

WTVALD         IS    A    FLAG    WHICH  IF    SF T    TO    1    CAUSES  VALIDATION 
INFORMATION    TO    BE   PRINTED    IN    ADDITION    TQ  THE 

LSUAL    OUTPUT    AT  EACH    PRINT    STEP.  OTHERWISE. 

NO    EXTRA    OUTPUT  IS  GENERATED. 


DECK  STRUCTURE. 


COMME  NT 

INDEX  • 
1 

2 
3 
4 


WHAT  FOLLOWS  FACH  INDEX  NUMBER  BELOW  MUST  START  A  NEW 
CARO. 


CARDS . 
GROUP  1. 
GRCUP  2. 
GROUP  3. 
SENTINEL. 

STL OOP    (FIRST  USE). 


FORMAT. 
SEE  A30VE. 
SEE  ABOVE. 
SEF  ABOVE. 
12 


STLOGP         IS    AN    INTEGER    VARIABLE    SET    TO  A 

SENTINFL    VALUE.     IF    STLOOP    IS    1    THEN  THE 
CARDS    THAT    FOLLOW    ARE    OBTAINED  BY 
RFCURSIVFLY    STARTING    AT    INDEX  3. 
OTHERWISE,    THE    NEXT    SFNTINEL    CARD  IS 
READ. 


SENT  INEL  • 

STLCOP    (SFCONC  USE). 


12 


STLOOP         IS    AN    INTEGER    VARIABLE    SF T    TO  A 

SENTINEL    VALUE.     IF   STLOOP    IS    1    THEN  THE 
CARDS    THAT  FOLLOU    ARE    OBTAINED  BY 
RECURSIVELY    STARTING   AT    INDEX  2. 
OTHERWISE,    THE    NEXT    SENTINEL    CARD  IS 
READ. 


SENTINEL. 

STLCOP    (THI RD   USE) . 


12 


STLOOP         IS    AN    INTEGER    VARIABLE    SET    TO  A 

SENTINFL    VALUE.     IF   STLOOP    IS    1    THEN  THE 
CARDS    THAT  FOLLOW    ARE    OBTAINED  BY 
RFCUSRIVELY    STARTING    AT    INDEX  1« 
OTHERWISE,    THE    RUN    IS  COMPLETED. 


EXAMPLE. 


35 


34  « 

C 

349 

C 

GROUP 

I    CARO( S) . 

350 

C 

GROUP 

2  CAPO, 

351 

c 

GROUP 

3  CARDS. 

352 

c 

1 

353 

c 

GROUP 

J  CARDS. 

354 

c 

0 

355 

c 

1 

356 

c 

GROUP 

2  CARD. 

ss-' 

c 

GROUP 

3  CARDS. 

358 

c 

>  0 

359 

c 

0 

360 

c 

1 

36  1 

c 

GROUP 

1    CARD! S) . 

36  2 

c 

GROUP 

2  CARD. 

363 

c 

GROUP 

3  CARDS. 

364 

c 

• 

36  5 

c 

• 

366 

c 

• 

36"' 

c 

0 

368 

c 

0 

369 

c 

0 

370 

c 

371 

c 

STRUCTURF 

OF    USER   SUPPLIED  SUBPROGRAMS. 

272 

c 

37  3 

c 

THE    DEFINITION   OF    THE   CALLING   PAHAMETFWS    USFO    IN    USER  SUPPLIED 

374 

c 

SUBRCUTINES   TOHV»CF    AND    DRAG   CAN   BE   PHTAINEO   FROM    THE  CURRENTLY 

375 

c 

SUPPLIED    VERSIONS   OF    THESE  ROUTINES. 

376 

c 

377 

378 

c 

379 

c 

SUBROUTINE   T  ORV*  OF  ( RPT  F ,    COSPTF,    SINPTF,    ZPTFt    VRWTFt  VAVnTF, 

380 

c 

•  ♦ 

VZWTF  «    WH I CH6 ) 

38  1 

c 

382 

c 

DECLARA  TI  CNS. 

383 

c 

384 

c 

GO 

TO    (100  .    500.    1000).  VKHICHG 

385 

c 

386 

c 

too 

SECTION   FOR   HEADING    IN    PARAMETERS    TO   BE   USF.D    IN  THIS 

387 

c 

SUBROUTI NE. 

3P8 

c 

389 

c 

KFTURN  ' 

290 

c 

39  1 

c 

500 

SECTION   FOR   COMPUTING    THE    VELOCITY    COMPONENTS   OF  THE 

392 

c 

TORNADO   WIND  FIELD. 

393 

c 

394 

c 

RETURN 

395 

c 

396 

c 

1000 

SECTION  FOR  PRINTING   ANY    RELEVANT    PARAMETERS   USED    IN  THF 

397 

c 

CCMPUT AT  ION. 

39  8 

c 

39  9 

c 

RETURN 

400 

c 

END 

401 

c 

402 

402 

c 

404 

c 

SUBROUTINE    ORAG{PARMTSt    DRAGCF  ,  wHICHGI 

405 

c 

36 


406 

C 

DECLARAT ICNS . 

407 

C 

408 

c 

GU    TO    (100,    500,     1000),  WHICHG 

409 

c 

410 

c 

100           SECTICN   f^CP    HEADING    IN    PARAMFTEHS    TO    BE   USED    IN  THIS 

41  1 

c 

SUtJRCUTI  NE. 

412 

c 

413 

c 

RETURN 

414 

c 

415 

c 

500           SECTION  roR  COMPUTING    THE    DRAG  COEFFICIENT. 

416 

c 

417 

c 

RETURN 

418 

c 

41<9 

c 

1000           SECTICN   FOR    PRINTING   ANY    RELEVANT    PARAMETERS   USED    IN  THE 

420 

c 

COMPUTATI  ON. 

421 

c 

422 

c 

RETURN 

421 

c 

END 

424 

c 

425 

426 

c 

427 

c 

PRCGRAf   USAGE  (OUTPUT). 

42« 

c 

42C} 

c 

LAYCUT    OF  OUTPUT, 

4  30 

c 

431 

c 

THE    PROGRAM   OUTPUT    IS    BROKEN    INTO    4  SECTIONS. 

432 

c 

433 

c 

SECTICN    1.      PROBLEM  DESCRIPTION, 

434 

c 

435 

c 

ThE    FCLLOWING    DESCRIPTIVE    INFORMATION    IS  PRINTED. 

43e 

c 

437 

c 

WIND   VELOCITY  PARAMETERS. 

43e 

c 

43q 

c 

ALL    THE    PARAMETERS    WHICH   ARE    INPUTEO   FOR    COMPUTED  THE 

440 

c 

TORNADO    WIND  FIELD    ARE    PRINTED.    THE    NUMBER,     NAMES,  AND 

441 

c 

MEANING    OF    THE    PARAMETERS    ARE    DETERMINED    BY    THE  USFR. 

442 

c 

443 

c 

DRAG   COEFFICIENT  PARAMETERS. 

444 

c 

445 

c 

ALL    THE    PARAMETERS    WHICH   ARE    INPUTED   FOR    COMPUTING  THF 

446 

c 

DRAG   COEFFICIENT    ARE    PRINTED,    THE    NUMBFR,    NAMES,  AND 

447 

c 

MEANING    OF    T(-E   PARAMETERS    ARE    DETERMINED    BY    THE  USER. 

448 

c 

44q 

c 

PARTICLE  PARAMETERS. 

450 

c 

451 

c 

SELf-F XPLANATORY .    FROM  INPUT. 

452 

c 

453 

c 

INITIAL   CONDITIONS    (TORNADO   WIND  FIELD). 

454 

c 

455 

c 

SELF-EXPLANATORY.    FROM  INPUT. 

45  6 

c 

457 

c 

INITIAL    CONDITIONS  (PARTICLE). 

45a 

c 

45<; 

c 

SELF-EXPLANATORY.    FROM  INPUT. 

460 

c 

461 

c 

SECTION    2.      TABULAR  RESULTS. 

46? 

c 

46  3 

c 

FIRST.    THE    INITIAL   TIME    AND   THE    TIME    INTERVAL    FOR    OUTPUT  ARE 

37 


4^  A 

Q 

PRINThO.    NFXT,     13    COLUMNS    OF    OUTPUT    AHF    GFNErtATtO.    EACH   L I NF 

465 

c 

CORRESPONDS    TO    A    SINGLE    PPINT    STEP.    THF    MEANING    OF  EACH 

t_  c 

c 

COLUMN  IS 

AS  FCLLOWS. 

46  7 

46  8 

c 

N  S  TP 

IS    THE    INTEGRATION    STEP    NUMBER   REACHED    AT  TIME 

4  ^  S 

T  . 

4  7 

c 

4  71 

T 

IS   THE   SIMULATED    TIME    AT    WHICH  OUTPUT  IS 

4  72 

c 

GENERA  TED • 

4  7  3 

4  74 

"■  XP 

ARE   THE    X,    Y,    Z    CCCRDlNATFS   OF    THE    PARTICAL  IN 

4  75 

YP 

THE    ABSOLUTE    FRAME    AT    TIME  T. 

4  7  6 

ZP 

4  7  7 

478 

c 

R  (  X  Y  ) 

IS    SQRT(XP    **    2    ♦    YP    **  2). 

4  79 

4R0 

c 

VXP 

ARE    THE    X.    Y.    Z   CCI^PQNENTS   OF    THE  PARTICLE 

48  I 

c 

V  YP 

VELOCITY    IN    THE    ARSOLUTE    FRAME    AT    TIME  T, 

4  B2 

c 

V  ZP 

4  P  ? 

4B  ^ 

HSP  EEO 

IS    SOHlT(VXP    ♦*    ?    ♦    VYP    **  2) 

4  3  5 

c 

4  R6 

SPEED 

IS    SqRT(VXP    **?■*■    VYP   ♦*    2   +    V7P   **  2) 

4  8  7 

4  R8 

c 

H  W  S  P  E  F  D 

IS    THE    HORIZONTAL    SPFFD    OF    THE    TORNADO   WIND  AT 

4r)9 

c 

THE    PARTICLE    IN    THF    ABSOLUTE  FRAME. 

490 

c 

49  1 

c 

RTF<XY ) 

IS    THF    HORIZONTAL    DISTANCE   OF    THE   PARTICLE  FROM 

492 

c 

THE    ORIGIN   OF    THE    TRANSLATING  FRAME. 

49  3 

c 

494 

c 

IF  WTVALO 

=    1    THEN    ADDITIONAL    OUTPUT     IS    GENERATED    AT  EACH 

495 

c 

PRINT  STEP 

.    SEE   THE    VALIDATICN   SECTION  FOR   THE   MEANING  OF 

496 

THIS   ADDITIONAL  OUTPUT. 

497 

c 

49  S 

c 

SECTION    3.      PROBLEM  TERMINATION. 

499 

c 

500 

c 

FIRST,  THE 

REASON    FOR    TERMINATION    IS    PRINTED.    IT    CAN    BE  FOR 

50  1 

Q 

ONE   OF    THREE  REASONS. 

E  02 

Q 

£  03 

Q 

1  )  THE 

FINAL    TIME    HAS    DFFK  REACHED. 

50  4 

Q 

505 

Q 

2  )  THE 

PARTICLE    HIT    THE  GRCUND, 

506 

c 

50  7 

c 

3)  THE 

ODE    SOLVER    FAILED  PPEMATURFLV, 

5  0  fi 

c 

5  0  9 

NEXT,    RESULTS    RELEVANT    TO    THE    *HCLF    RUN   ARE    PRINTED.    EACH  IS 

510 

c 

CURRENTLY 

TAKEN    TO    BE    THE    MAX  IMUM    IN    ABSOLUTE    VALUE  OVER 

511 

c 

VALUES  GENERATED   AT    THE    PRINT    STEPS.    THEY    HAVE  THE  FOLLOWING 

512 

MFANI  NG. 

513 

514 

c 

MAXVXP 

ARE   THE   MAXIMUM    X,    Y,    7    COMPONENTS    IN  ABSOLUTE 

515 

c 

MAXVYP 

VALUE    OF    THE    PARTICLE    VELOCITY    IN   THE  ABSOLUTE 

516 

MAXVZP 

FRAME  . 

5  17 

518 

c 

MA  XMVP 

IS    THE    MAXIMUM    SPEED    OF    THE    PARTICLE    IN  THF 

51  9 

c 

ABSOLUTE  FRAME. 

520 

c 

52  1 

c 

SFCTIGN    4.  PLOTS. 

38 


J 


! 

i 


f- 

=  0  T 

c 

/FRfl    DR    MPRF    PI  D  T AHt-     PRnDlirFP    f  1 NF     PFR    PARF«     T  HP  HRDlNATt 

524 

AND    AdSCISSA    LAHELS    CAN   gF    FCUND    CENTFPED    BFrLOW    THE  GRAPH. 

525 

c 

THE:    PLOTTING    CHARACTF.R    IS    THE    LETTER    X.    AS    NOTED    IN  THE 

526 

SECTION.    APPLICABILITY    AND    RE  S  Tk  I  C  T  10  N  S  ♦    AT    MOST    101  DATA 

527 

c 

POINTS    ARE  PLCTTED. 

528 

c 

c  o  c 

f 
\^ 

«  O  U 

c 

1 

7  O  I 

f 

INPUT  rJ-FTKINR. 

53  3 

THE    FOLLOWING    VARIABLES    ARE    CHECKED    FOt^    VALIDITY    IN    THE  MAIN 

534 

PROGRAM*     IF    AT    LEAST    ONE     IS    INVALID    AN    ERROR    MESSAGE    IS  PRINTED 

535 

AND    THFY    ARE    ALL   RE SF T    TO   DEFAULT    VALUES.    THE    DEFAULT  VALUES 

536 

CAS   BF   FOUND    IN  THE    SECTiCN.    MACHISE/SYS TFM    DEPENDENT  FEATURES* 

£3T 

c 

AT    THE    START   OF    THE    EXECUTABLE    C  TDE    IN    THIS    MAIN  PROGRAN'. 

53  8 

539 

HI                   IS    THE    INITIAL    STEP.    IT    MUST    BE    POSITIVE    AND  LESS 

540 

c 

THAN   FXTM I T. 

541 

542 

Q 

FXTMIT         IS    THE    TIME    INTERVAL   USFO   FOB    OUTPUT.    IT    MUST  BE 

543 

pns  IT  I  VF  a 

£4  4 

£45 

c 

PROGRAM    DETECTABLE  ERRORS* 

£46 

547 

c 

AFTFR   EACH    INTEGRATICN   STEP   THE   VARIABLE    lER    IS    CHECKED.  IF 

54B 

c 

ITS    VALUE    INDICATES    THAT    THE    LAST    STEP   FAILED    THEN  INTEGRATION 

549 

c 

TERMINATES    PREMATURELY*    THE   REST   OF   ThF    OUTPUT    IS  GENERATED 

b:<%  n 

DO  u 

55 1 

c 

WILL    REFLECT   ONLY   RESULTS    GENERATED    UP    TC    THE    TIME    OF  THE 

5  52 

P  C  ROD  . 

5  5  3 

554 

DisrussiuN  nF  method  and  ai  gorithm. 

ec;5 

556 

c 

THE    PROBI  EM    Tn    HE     SnLVFD    IS    A    NUMFPICAL     INITIAL    VALUE    PRTHLEM  IN 

557 

ORDINARY    DIFFERFNTlAl     EQLjATION"^.     SIX.    FTR^T     CHDFK.  ORDINARY 

556 

<- 

DIFFFRFmTIAI       FOUAT  inN^.     HF  PRF  SFN  T  I  Nf;     Th  f-     FOUATfnN^    flF     MflTfriNl.  ARF 

559 

c 

INTEGRATED    FROM    SOMF     I  N  I  T  I  AI     T  1  MF    WITH    SPFCIFTED  INITTIAI 

ft  •  ^  %  1  >«J  r»        1  u_           i                            U  >^U—        ft  I  ^  ft    I    ft  ^  ^_       1    ft        i_      Vf   ft  |   •  t       j       i  _  \^  ft  ■     i  i_  (_/       ft  TN  ft    •    ft  ft~ 

560 

Q 

misjniTIONS    TO    SnMF    LATFR    time     V^HICH    SATISFIES    SOMF  TFRWINATinN 

561 

rCNO I T I C  N. 

Va  ^   1  ^             ft       1       ft    U    1  ^  V 

562 

563 

c 

THE    ODE    SOLVER*    VOADAM*     IS   BASFD    ON    THE    VARIABLE    ORDER*  VARIABLE 

564 

STEP    ADAMS    fPCF    CFSirNFT    ANn     I MPI  EMFNTATFD    PY    C.    W.    GFAR  IN 

565 

RFFFRFNCF     1  )  .     S I  Nf"  F     THF-     tiV'^TFM    OF    ODFS    TP    HF    «ini  VFD     I        NFN  — ^TIFF 

566 

AI   1      PARAMFTFRS    AND    r"  AI  1  S    TO    SURRHUTINFS    RFOI   IRFO    FOR     c;ni   U I  ^slr; 

567 

STIFF     SYCTFM^    HAUF    RFFN    OFl  FTFO. 

568 

/- 

569 

CN    EACH    CALL    TO    VOADAM    THF    DDF    SOI  VFR     IS    ASKED    TO     INJTFGRATF  THF 

570 

^YSTFM    of    onFS    nVFR    a     STFP    Of    I  FNiHTt-i    H.     THF     V  AI  lIF    OF     H    AWI'iF"^  FROM 

571 

CNF    flF    T  l«n  SOllRrFS. 

572 

573 

Q 

1  )     THF    VAI  liF     nF    H    RFTIIRNFD    RV     THF    PRFv/IHl  S    C  AI  1      OF  VOADAM. 

574 

575 

c 

?  )     TH  F    VA 1  (JE    CF    H    SPFCIFIFD    Hv    CAIIFR     T  N  TFR  ACT  I O  M  . 

£  9       1  n          w  ^  w  w  t.      ^«        'I     j~tv»ftFfti_L'     '.JT  *— ■     ^  1  t_          i  IN  I  L      ^  ^  1    ^  u  IN  • 

576 

577 

c 

SOURCE    2)    IS    USED    ONLY   FOR    THF    INITIAL    STEP    AND    WHEN    THE  STEP 

578 

c 

IS    MODIFIED   SC   THAT    IT    FALLS    ON    A    MULTIPLE    OF    THE    FIXED  TIME 

579 

c 

INTERVAL   USED  FCR    PRINTING.    CTHEfiwISE,    SOURCE    I)    IS   USED   SO   AS  TO 

39 


50  0 

C 

ACHIEVE    AN    FCCNCMICAL     I  NT  E  G^^  AT  lO  N  . 

56  I 
58  ? 

C 

c 

THE  INPUT 

VALUE    OF    H    IS    USED    UNLESS    THE    FQPCR   CRITERIA    CANNOT  RF 

5  H  3 

c 

MET.    IN    THIS    CASE    THE    STEP    AND/ CP    OfiDEW    AWE    MDCIFIED    TO    TRY  TO 

5  84 

c 

MEET  THt 

ERROR    CRITERIA,    IF    AN    ATTEMPT    IS    MADE    TO   REDUCE    THE  STFP 

585 

c 

BELOW    A    CALLER   SUPPt.  IFO   VALUEt    HMIN,    rnf    (JOE    SOLVER    QUITS  AND 

5P6 

c 

RETURNS    THE    APPROPRIATE    N0N7ER0    EHKCW    CINOITION  CODE. 

58  7 

c 

58  8 

c 

ONCE    A  SUCCESSFUL    STEP    HAS    HEEN    TAKEN,    VOAOAM   ESTIMATES  AND 

5fly 

c 

RETURNS  A 

GCOD    VALUE    OF    h    TO    BE    USED    FOR    THE    NEXT    STEP.  THIS 

590 

c 

EST IMATED 

STEP    CANNOT   RE   GREATER    THAN    A    CALLER    SUPPLIED  VAtUF, 

591 

c 

HMAX  . 

59;? 

c 

59  3 

c 

FOR  MORE 

DETAILS   UN  VOAOAM  CONSULT    THE   MACHINE  READABLE 

5  94 

c 

OOCUMLNTATI  GN    AT   THE    BEGINNING   OF  THE  SUBROUTINE. 

595 

c 

£  9  6 

c 

VAL lUAT ICN. 

59  7 

c 

59  fi 

c 

THE    PWOGRAM    FRCVIDES   THE   USER    BY    WAY    OF    THE    WTVALO    INPUT  VARIABLE 

599 

c 

THE    CAPABILITY    OF    PRINTING    IMPORTANT    INTERMEDIATE  QUANTITIES. 

f>0  0 

c 

THESE    QUANTITIES    ARF   PRINTED    AS    ADDITIONAL    OUTPUT    AT    THE  NORMAL 

6  01 

c 

PRI  NTI  NG 

STEPS.    THE    AMOUNT    OF    OUTPUT    PER    PRINT    STEP    INCREASES  FROM 

60  ?. 

c 

ONE  LINE 

TO   SEVEN   LINES.   TO    INTERPRET    THE   MEANING  OF    THE  VARIABLES 

eo3 

c 

PRINTED   CONSULT    THE   TABLE    BELOW    ACCOMPANIED    BY    APPENDIX    A  OF 

604 

c 

REFERENC  E 

2  >  . 

6  0S 

c 

606 

c 

R  TF 

IS    THE    CYLINDRICAL    RADIUS   OF    THE   ORIGIN   OF  THE 

607 

c 
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TORNAOC*W  INC (  1  )  .TCKWOF  $FX ( 0  ) 

SUBRCU  T  I  NF    T C K WOF ( fiP T F , C 03 PT  F , S INPTF . ^PTr,V«WTF,VAWTF,VZ*TF. 
»  WHICHO) 

THIS    SLBRCUTINE    COMPUTES    X  t-f:    TCRNADO    V.tNO   VtLOCITY  COMPONFNTS 

DFFIMTICNJ   OF  PARAMETERS. 

RPTF .COSPTF, SI NPTF , ZPTF         =    RADIAL    CCCRDINATB,    COSINF    AND    SINE  OF 

THE    ANGLL  ,    AND    Z    CCCROINATE    OF  THE 
^ARTICLE    WITH   RESPECT    TO    THE  CYLINDRICAL 
COORDINATE    SYSTEM    DEFINED    AS  ThF 
J      .  TRANSLATING  FRAME 

VR WTF , VA WTF , VZWTF  =    RADIAL,    ANGULAR,     AND    Z    COMPONENTS   OF  THE 

WIND    AT    THE    PARTICLE    AS    MEASURED    IN  THE 
TRANSLATING  FRAME 
WHICHG  =    1     -    READ    IN    WIND    VFLOCITY  PARAMETERS 

■  .  ,  2    -   COMPUTE    '*IND    VELOCITY  COMPONENTS 

3    -   PRINT    WIND   VELOCITY  PARAMETERS 

PARAMETER  DECLARATIONS. 

REAL    RPTF,  COSPTF.  S  INPTF  ,ZPTF  ,  V  R  W  TF  ,  V  A  W  TF  ,  VZ  Vn  TF 
INTEGER  \»HICHG 

CTHER    COMMON  VARIAGLES. 

REAL    RMTVZO .MTV.K  I  ,K?,K 3.K4, K5 

COMMON    / TWDFDR/RMT VZO .MT V ,K 1  ,K2  ,  K3  ,  K4 , K5 

LOCAL    VARIABLE  DECLARATIONS. 

REAL    RMTVZ  ,R33Z  ,RAT  10 
INTEGER  ID.QO 


MACHINE/SYSTEM    DEPENDENT  FEATURES. 

DEFINITION    OF    I/O    UNITS    USED    IN    THIS  SUBROUTINE. 
10  IS    THE    UNIT    USED    FOR    INPUTING  DATA. 

CD  IS    ThE    UNIT    USED    FOR    OUTPtJTiNG  RESULTS. 

DATA    ID.    OD/5  ,  6/ 

GO    TO    (  100,500,  1000>  ,«^H1CHG 

READ    IN    WIND    VELOCITY  PARAMETERS 

100  REAC(IC.250)    RMTVZO  ,  MTV  ,K  1  ,  K2  ,K3  ,K<»  ,  K5 

250  FORMAT( 1 CES. 0) 

RETURN 


I 

2 

3 

C 

4 

C 

C 

6 

C 

7 

C 

8 

C 

C 

1  0 

C 

1  1 

c 

1  2 

c 

1  3 

c 

1  A 

c 

1  5 

c 

1  6 

c 

1  7 

c 

1  e 

c 

1 

c 

20 

c 

2  1 

2  2 

23 

c 

24 

c 

25 

c 

26 

2  7 

c 

28 

2<5 

c 

30 

c 

3  1 

c 

32 

34 

c 

35 

c 

36 

c 

3  7 

c 

38 

c 

3<; 

c 

40 

c 

4  1 

c 

42 

c 

43 

c 

44 

c 

45 

c 

4  e 

47 

c 

4  8 

c 

4S 

c 

50 

c 

51 

52 

c 

5  3 

c 

54 

c 

55 

56 

57 

42 


56 
5Q 
60 
6  1 
f>2 
63 
64 
^  5 
66 
67 
6°l 
69 
70 
71 
72 
73 
74 
75 
76 

78 
79 

80 

e  t 

82 
B3 
84 
85 
86 
87 
88 
89 
90 
91 


100 
1  0  1 
102 
END  PBT 


COMPUTE    WIND    VELOCITY  COMPONENTS 


500 


525 
550 


575 


IF( ZPTF . GE . 60 . 0 )    GO    TO  525 
RMTVZ=RWTVZO ♦Kl ♦ZPTF 
GC    TO  550 

RMTVZ=RMT VZO+K I ♦60.0 

P33Z=  (  WTV/33  .0  )  ♦♦  .625*RMTVZ 

IF (RPTF. GE  .RMTVZ »    GO    TO  575 

NATI 0=RPTF/RMTVZ 

GO    TO  6C0 
HAT IC  =  RMTV2/RPT  F 
IF ( RPTF.GE.R33Z )    GO    TO  650 


600  IF  (  ZPTF  ,  GE  .60  .0  >    GO    TC  625 

VR  WTF=-K5^(R33Z-RPTF )/(HJ3Z- RMTVZ) *RPTF 
VAV»TF=R  AT  IO*MT  V 

V7V»TF=K3+  {R33Z-RPTF)/  (  R3  3Z-PMTVZ  )  ♦ZPTF  ♦K 4/ 3  .0  ♦VAWTF 

RETURN 

625  IF( ZPTF. GE. 240.0)    GO   TO  650 

VRHiTF=-K5^(  R33Z-RPTF  )  /  (R3  3Z-KMTVZ  )  ♦PP  TF  ♦(  240.  0-ZPTF)  /  I  80. 
VA*TF=RATI C+WTV 

VZVgTF=(  l,3  3-ZPTF/180.0)^(K3^(F!J3Z-NPTF)/(H3  3Z-RMTVZ)^60.0 
*  K4/3  .0  ♦V  AWTF  ) 

RETURN 
650  VRi«TF=0.0 

V  AV>TF  =  RAT  lO^MTV 
VZ  ViiTF  =  0.  0 
RETURN 

PRINT    WIND   VELOCITY  PARAMETERS 


1000 


VKRI  TF  (  00  .  1  1  00  )    WMTVZO  ,MTV  .Kl  .K2,K3.K4,K5 


92 

1100 

FORMAT! 26H0W INO 

VELUC ITY 

PARAMF  TE 

93 

* 

6X.15H  RMTVZO 

.  1PE12.4/ 

94 

♦ 

5X.15H  f'TV 

,  IPE I  2. 4/ 

95 

♦ 

5X  .  1 5H  Kl 

. IPE 12. 4/ 

96 

* 

5X.15H  K2 

, IPEl 2.4/ 

97 

♦ 

5X,  1 5H   K  3 

. IPE 12.4/ 

98 

* 

5X . 1 5H  K4 

,  IPEl 2  .4/ 

99 

* 

£  X  ,  1  5H    K  5 

.IPE12.4  ) 

RETURN 


FND 


fflPRT ,5  F 1 ,CRA6*EX 
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TOPNAOn*  Vk  I  NU  (  1  )  .  OPAG  Eh  X  (  0  ) 


1 

oUnKUU  T  I  Nr     L>RAC>fHAKM  lit  UN    t>v  r  •MMlCMol 

C 

3 

c 

* 

c 

c 

c 

AiTCC-'^T     T  Lj  1-      /"OAJDllT  AT  ItfltJ     ntT     T  MP      HO  A 

P 

c 

t*lrrr  It-  IrNT 
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I  0 

I  1 

uiM  1  ri-4                                                             =     1     -    RF  An     T  M    DP  Afl    rnf  i-  FiriKNiT    PARA  MF  TL 

c 

?   —  rnMPiJTF  DRAH  rnfFFiriFwT 

£              \.  X  J  ■'1     \j  1  r              r4          \«'Jl«»        IV-       -  • 

I  3 

i. 

"I    —    PRIMT    P'RAT*    rnFPFf/'IFMT    PA  OAMFTPK^ 

1  A 

I  *» 

f 

1  ~ 

c 

DADAMFTFP     HTPl  ADATrrK.-^. 

1  A 

I  f> 

V. 

1  f 

LJCAlCACJfciTC/l% 
HtAL     HAKprlT^ll  f 

I  fi 
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I  V 

INit-tjrHwr'iCrivj 

o  t 
£  1 

V- 

fiTi-iFD    ^riMiMnKj    WAtJiAm  f^q^ 

C.  w 

DF  A  1      r  I"^© 

O  A 
c  *♦ 

^  C 

Km 

O  T 

Km 

1   P^AI       WADfAPl    F     PFP*!  AOATIflUQ^ 

TKiTPPF-D      ir*  no 

3u 

K. 

J  1 

K. 

f 
V. 

UA/^LjfKtCyCVCTCU      nCOCTklOClklT  STCATllLiPC 

C 

Ob 

V. 

r^PiriKi  TfOKj             fyp    iikjitc    iicPn    ik'    TUtc  i^iic^uniiTffAiP 
Ut  PllNlllUiN    Up     l/U    UNIlD    CIjpLI     ll\     irili>    oUOKUlJ  T  1 NL.  • 

K. 

3  7 

c 

lU                      lb    T  Ht.    UNIT    UStiU    rUK     INPUTING    O  A  1  A  • 

3  e 

39 

K. 

L  U                         I  3     1  Ml:.     UN  IT     U  oP  L'    rUK     U^^TKU'l  f^t>     Kr.  OUL' 

4  C 

c 

4  1 

c 

4  ^ 

□  ATA     lC»     OD/^TS*  O/^ 

4  3 

c 

4  4 

C 

A 

A  o 

C- 

A  £ 

4  C 

I. 

4  r 

TU     ilUOtbUUf  iCUUlt  AM  1  CM(j 

A  a 
4  tJ 

\. 

A  O 

K  t.  AU     IN    UH  AO    C-Utr  p  It.  IrNT  HAKAMcTpHS 

&0 

c 

5  I 

\  />  f\              oCAPk/if^  ^r\f5A/™ 
luO             WtAUilUt250I  t-UHAvj 

52 

2  O  CJ            p  UHMA  T  I  or.  1  ^  •  O  1 

53 

54 

c 

55 

c 

COMPUTE    0«AG    COEFFIC lENT 

56 

c 

57 

50  0           DRA  GCF-CDRAG 
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5P 

RFTLRN 

c  q 

c 

60 

c 

PRINT    DPAG   COEFFICIFnT  PARAMETERS 

t  I 

c 

6? 

1000           «tRITE(OC«l  100  )  CORAG 

6  3 

1100           FORMA  T(  29H0DRAG    COEf-FICIENT  PARAytTEHS./ 

64 

*                   5X,15H   CDRAG                   =  ,1PE12.A> 

e  = 

RETURN 

66 

c 

67 

END 

END  PRT 

aPPT.S  Fl.CATASEX 
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TCRNADQ*teIKO(l)  .CA>TA*fcX(l  ) 


1  46. C         130.0  0.0  0.0  1.0  1.0  1.0 

2  1  .0 

3  1.0  10.0 

4  46.0  0.0  40.0  0*0  0.0  0  .0 

5  0.0  0.0  35.0  0.0  0.0 

6  10.0 

7  l.Or-4  l.OE-1  l.OE-4 

a  1  1  1  1  1  I  0 

s  c 

10  c 

11  0 


END  FRT 
a)XQT    F  I  .ABS«E  X 
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